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FOREWORD 


NASTRAN® (NASA STRUCTURAL ANALYSIS) is a large, comprehensive, 
nonproprietary, general purpose finite element computer code for structural 
analysis which was developed under NASA sponsorship and ' became ^ailable to 
the public in late 1970. It can be obtained through COSMIC® (Computer 
Software Management and Information Center), Athens, Georgia, and is wi y 
used by NASA, other government agencies, and industry. 

NASA currently provides continuing maintenance of NASTRAN through COSMIC. 
Because of the widespread interest in NASTRAN, and finite element ^ 1 ?°* the 
general, the Twentieth NASTRAN Users' Colloquium was organized and held at the 
Sheraton Colorado Springs, Colorado Springs, Colorado . on April 27 - 
1QQ? (PaDers from previous colloquia held in 1971, 19/*:, *-9/o, » ’ 

1977’, 1978^ 1979, 1980, 1982, 1983 1984 1985, 1986, 19 f' 2893 

and 1991 are published in NASA Technical Memorandums X-2378 . X-2637. X 2893, 
X-3278 X-3428* and NASA Conference Publications 2018, 2062, 21di, cl y, 

2284 2328 2373, 2419, 2481, 2505, 3029, 3069 and 3111.) The Twentieth 

Colloquium ’provides some comprehensive general papers ° n h aoorolches 

finite element methods in engineering, comparisons with other a PP r ° a ches, 
unique applications, pre- and post-processing or auxiliary programs, and new 
methods of analysis with NASTRAN. 

Individuals actively engaged in the use of finite elements or NASTRAN 
were invited to prepare papers for presentation at the Colloquium. ese 
naoers are included in this volume. No editorial review was provided by NASA 
or^COSMIC; however, detailed instructions were provided each author to achieve 
reasonably consistent paper format and content. The opinions and data 
presented 7 are the soleresponsibil ity of the authors and their respective 
organizations . 


NASTRAN® and COSMIC® are registered trademarks of the National Aeronautics and 
Space Administration. 


iii 


PRECEDING PAGE BLANK NOT FILMED 



CONTENTS 


Page 


FOREWORD . 


1. NASTRAN INTERNAL IMPROVEMENTS FOR 92 RELEASE 

by Gordon C. Chan 

(UNISYS Corporation) 

2. EVOLUTION OF A NASTRAN TRAINER 

by H. R. Grooms, P. 0. Hinz and M. A. Collier 

(Rockwell International) 

3. ANIMATION OF FINITE ELEMENT MODELS AND RESULTS 

by Robert R. Lipman 

(David Taylor Research Center) 

4. ACCURACY OF THE TRIA3 THICK SHELL ELEMENT . . . • 

by William R. Case, Marco Concha and Mark McGinnis 

(NASA/Goddard Space Flight Center) 

5. VALIDATION OF THE CQUAD4 ELEMENT FOR VIBRATION AND SHOCK ANALYSIS 

OF THIN LAMINATED COMPOSITE PLATE STRUCTURE 

by Douglas E. Lesar 

(Naval Surface Warfare Center) 

6. A CASE OF POOR SUBSTRUCTURE DIAGNOSTICS 

by Thomas G. Butler 

(Butler Analyses) 


7. ALTERNATIVE METHODS TO MODEL FRICTIONAL CONTACT SURFACES USING 

NASTRAN 

by Joseph Hoang 

(GE Government Services) 

8. HIERARCHICAL TAPERED BAR ELEMENTS UNDERGOING AXIAL DEFORMATION 
by N. Ganesan and S. K. Thamoi 

(GE Government Services) 

9. TRANSIENT THERMAL STRESS RECOVERY FOR STRUCTURAL MODELS • • • • 

by William Walls . 

(McDonnell Douglas Space Systems Co.) 

10 TRANSIENT LOADS ANALYSIS FOR SPACE FLIGHT APPLICATIONS . . . . 

by S. K. Thampi, N. S. Vidyasagar and N. Ganesan 
(GE Government Services) 


116 

124 

134 

154 


11. ACOUSTIC INTENSITY CALCULATIONS FOR AXISYMMETRICALLY MODELED 

FLUID REGIONS •••;••* 

by Stephen A. Hambric and Gordon C. Everstine 
(David Taylor Research Center) 


v 


pFfECEOING PAGE 


blank NOT FILMED 




N9 2- 24 325 

NASTRAN INTERNAL IMPROVEMENTS FOR 92 RELEASE 


by 

Gordon C. Chan 
NASTRAN Maintenance Group 
UNISYS Corporation 
Huntsville, Alabama 


SUMMARY 

The 1992 NASTRAN release incorporates a number of improvements transparent to users. 
The NASTRAN executable has been made smaller by 70 percent for the RISC base Unix 
machines by linking NASTRAN into a single program, freeing some 33 megabytes of system 
disc space that can be used by NASTRAN for solving larger problems. Some basic matrix 
operations, such as forward-backward substitution (FBS), multiply-add (MPYAD), matrix 
transpose, and fast eigensolution extraction routine (FEER), have been made more efficient by 
including new methods, new logic, new I/O techniques, and, in some cases, new subroutines. 
Some of the improvements provide ground work ready for system vectorization. These are finite 
element basic operations, and are used repeatedly in a finite element program such as 
NASTRAN. Any improvements on these basic operations can be translated into substantial cost 
and cpu time savings. This paper also discusses NASTRAN in various computer platforms. 


NASTRAN SINGLE LINK 

The 91 NASTRAN RSIC base Unix version was released as a single link program. (The 
multi-link version can also be built.) The shell program that controls NASTRAN s execution was 
modified so that it can run NASTRAN under single- or multi-link environments automatically. 
The success of the single link Unix version (previously referred to as NASTRAN superlink) 
prompts the conversion of the 92 VAX/VMS NASTRAN into a single link program. The VAX 
single link program is extremely useful in debugging the Unix RISC version since the latter has 
extremely poor system diagnostics and provides no error trace-back. The biggest advantage of 
the single link version, and particularly in a small workstation environment, is that it needs only 
25 percent of the disc space to hold the executable program as compared to the multi-link 
version, and 33 megabytes of disc space is returned to the system that can be used by NASTRAN 
for solving larger structures. It is also faster by almost 70 percent to build the single link 
NASTRAN as compared to the complete multi-link version. In program execution, the single link 
has no overhead for link-switching, and data saving and recovering between links are not needed. 
The single link versions, both in RISC base machines and VAX/VMS, appear to be running faster 
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since the program could be executing several "links" ahead of the screen printout if NASTRAN 
is executed interactively. For purposes of consistency, the single link program issues the link 
numbers as if the program were generated in multiple links. Presently, there are no single link 
versions for IBM, CDC, and UNIVAC machines. 


NASTRAN IN VARIOUS MACHINE PLATFORMS 

COSMIC/NASTRAN is supported in five computer platforms: IBM, CDC, UNIVAC, 
VAX/VMS, and DEC/ULTRIX (a RISC base Unix machine). However, the COSMIC/NASTRAN 
is continuously improved aiming towards a unified environment. ANSII standard FORTRAN 77 
is used; and most of the machine dependent items are removed if possible. Some source codes 
are modified or re-written with system vectorization in mind. A far-reaching plan is initialized 
in the boot-strap (BTSTRP) subroutine that is used to set up machine-dependent constants for the 
five machines that COSMIC supports. The 92 COSMIC/NASTRAN has expanded the 
machine-dependent constant table in BTSTRP for 18 major computers and two dummies. A user, 
or a third-party organization, can install NASTRAN to a new platform currently not supported 
by COSMIC. Also, the users or third parties using the new computer platform can talk to one 
another since the machine type has been pre-arranged. The machine-constant table in BTSTRP 
is arranged for the following computers: 

DUMMY, IBM, UNIVAC, CDC, VAX, DEC/ULTRIX (RISC base Unix), SUN, 

IBM/AIX, HP, SILIC.GRAPHICS, MAC, CRAY, CONVEX, NEC, FUJITSU, 

DATA GENERAL, AMDAHL, PRIME, 486, DUMMY 

BTSTRP sets up the machine constants correctly only for the first six machines. The first 
DUMMY is set up for IBM 7094, which has long been obsolete, and therefore can be reused for 
any machine not on the list. The last DUMMY is intended for the same purpose. The constants 
for the remaining 14 machines are dummies or guess- values. Therefore before moving 
COSMIC/NASTRAN to any computer platform, one must first re-supply the correct machine 
dependent constants in BTSTRP for that machine. The machine constants are well defined in the 
BTSTRP subroutine, such as NBPW, number of bits per woid; NBPC, number of bits per 
character, NLPP, number of lines per printout page, etc. 

There are also a few machine-dependent constants outside BTSTRP that need to be set 
locally. These few constants are scattered in about eight or nine NASTRAN subroutines. For 
example, some of these constants involve convergent criteria used only locally, and therefore not 
at the same level of importance as those in BTSTRP. To locate these few machine-dependent 
constants has been made easy in the 92 NASTRAN release. If the first word of labeled common 
/MACHIN/ in a subroutine is MACHX, not just MACH, there are machine dependent constants 
that require fixing. 

The main NASTRAN program in the single link environment is called NASTRN. In the 
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multi-link environment, each of the 15 links is a complete program by itself, and the main 
programs are NAST01, NAST02, NAST03,..., NAST15. There is something very special in 
NASTRN and in NAST01 in the 92 NASTRAN release. If the DEBUG flag in NASTRN or 
NAST01 is changed to 1, the so-called <LINK 1> portion of the single link NASTRAN, or the 
regular LINK 1 in the multi-link program, will go through a series of machine compatibility 
checks The results will echo back to the user. If some things, or some parameters, are set 
incorrectly, NASTRAN will stop. For example, if the FORTRAN OPEN statement of a new 
computer platform is in byte count for the record length, RECL, and the user sets the 
corresponding constant in BTSTRP to word count, an error diagnostic will appear. 


Due to COSMIC policy, the NASTRAN Maintenance Group has never tried to move 
NASTRAN outside the five designated computer domains. However, it was demonstrated by a 
user in 1991 that the DEC/ULTRIX version required only two minor changes to move 
NASTRAN to a SiliconGraphics machine; and those two changes have been incorporated back 
into the 1992 COSMIC/NASTRAN release. The NASTRAN Maintenance Group in UNISYS 
welcomes any contribution from the users or third parties, concerning the migration of 
NASTRAN to different computer platforms. 


INTERNAL IMPROVEMENTS 

The PROFILER, a performance analyzer of the VAX/VMS machine, was used to identify 
the major time consuming elements in several typical NASTRAN runs. This was followed by a 
major effort to look into the logic, computing mechanism, methods of calculation, order ot 
execution, paging, vectorization etc. of the time consuming areas, and to search for improve- 
ments. Most of the time consuming elements are basically standard matrix operations, and they 
are the essential elements of a finite element program. Their treatments are either textbook 
standard”, or "company proprietary". Generally speaking, the original NASTRAN developer did 
a very fine job in these areas. Further improvements are not easy, and cannot be treated lightly 
or as short projects. Indeed, several of the 1992 internal improvements were made in periods of 
several weeks, not days. Some improvements were made on top of previous improvements or 
newly written subroutines. Sometimes a simple line of improvement may require days of careftil 
study, thorough understanding of the program algorithms, and the theoretical treatment of the 
subject The final improvements in the 92 NASTRAN release are quite satisfactory. Several slow 
moving areas are speeded up 25 to 30 percent, and in some areas, three to five times faster 
Appendix A tabulates some of the test results. Most of the new internal improvements of the 92 
NASTRAN release can be removed by DIAG 41. 

The following sections of internal improvements apply only to large matrix operations 
Usually the matrices are many times larger than the available computer core memory can hold 
at one time. Many of the matrix operations must be done by parts, or in a number of passes. 
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IMPROVED FORWARD-BACKWARD SUBSTITUTION 


The forward-backward substitution (FBS) is used for matrix inversion, load-solution, 
eigenvalue iteration, and many other applications that follow matrix decomposition. If the number 
of columns on the solution side of the equation is large, FBS can be very time consuming. It 
could easily take 10 to 20 times longer to go through FBS than to do the matrix decomposition. 

In NASTRAN, the driver for FBS is the subroutine FBS. The actual FBS computation 
takes place in FBS1, FBS2, FBS3, or FBS4 for real single precision, real double precision, 
complex single precision, and complex double precision calculation respectively. If FBS requires 
more than seven passes, the new improvements will automatically kick in. The improvements are 
in four new subroutines, FBSI, FBSH, FBSHI, and FBSIV, similarly arranged as the subroutines 
FBS1/2/3/4. The new improvements include reduced I/O operations, large data blocks, and new 
row-and-column matrix multiplication. The new FBSI/H/IU/IV are about 30 to 50 percent faster 
than the original FBS1/2/3/4, as tested in COSMIC’s VAX/VMS machine. 


THE FEER METHOD 

Seventy to seventy-five percent of the cpu time used by the FEER method (the fast 
eigenvalue extraction method, with real tridiagonal reduction) is actually spent in the subroutine 
FNXTVC (double precision version) or FNXTV (single precision version). The main time 
consumer in FNXTVC/FNXTV is the forward-backward substitution operation. Untike the FBS 
module, the open core in FEER method is not fully used, and particularly not in FNXTVC/ 
FNXTV subroutines. The improvements in FNXTVC/FNXTV include reduced I/O operations, 
full utilization of the core space, and new row-and-column matrix multiplication. The new 
improvements in FNXTVC/FNXTV alone produce impressive results - reducing the FEER 
method cpu timing by 30 to 200 percent, as tested on the VAX/VMS machine, and on a CRAY 
(tested on RPK/NASTRAN). 

COSMIC/NASTRAN sometimes gives negative values for the rigid body eigenvalues. 
(They should be zeros.) Sometimes the negative values could be quite large (-1.E+5 range) 
particularly on the IBM machine. The explanation for this strange behavior, and the solution to 
the problem, may or may not match. On the solution side, since the rigid body frequencies are 
zeros, the 92 COSMIC/NASTRAN FEER method, by default, will set them to zeros. The second 
solution option is to reinforce certain key areas of computation in FNXTVC/FNXTV by 
quad-precision (real* 16) for the 32- bit word machines, and by double-precision for the 60- or 
64-bit word machines. This second option gives good results and moves the rigid body 
frequencies down to l.E-6 to l.E-12 range. However, it takes 2 to 3 times longer to compute. 

To activate quad-precision calculation in a 32-bit word machine or double-precision in 
a 60- or 64-bit word machine in FEER method, one only needs to replace the BCD word "FEER" 
on the third field of the EIGR bulk data card by "FEER-Q". Replacing "FEER" by "FEER-X" 
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will prohibit the substitution of zeros for the rigid body frequencies. 

Now to explain what happens to produce the negative eigenvalues. The criterion for 
orthogonality convergence EPSILON is l.E-28 for the 32-bit word machines using double 
precision computation, and l.E-24 for the 60- and 64-bit word machines using single precision 
computation. In FNXTVC/FNXTV, the accumulated sum of a mode shape is compared to 
EPSILON. The accumulated sum at the end of a do loop is the difference of two very close 
numbers. Therefore this very small numeric difference is a function of the number of digits a 
computer word can hold. Most 32-bit word machines with double precision calculation, and most 
60- or 64-bit word machines with single precision, are limited to 14 to 16 digits per word, and 
therefore down to only l.E-14 to l.E-16 numeric accuracies. Here, the mantissa as well as the 
exponent of a data word is important! Since all five computers supported by COSMIC exhibit 
this accuracy problem, it becomes meaningless to do an orthogonality check by comparing the 
result to EPSILON, which is another 10 to 14 decades smaller. To home in to the size of 
EPSILON could end up producing big numeric errors. (Just a guess here.) Since the mantissa of 
an IBM machine word, in double precision, is smaller than a VAX, IBM seems to produce big 
negative rigid body eigenvalues more often than the VAX. Similarly the 60-and 64-bit machine 
using single precision computation could be worse than IBM. 

In the 92 NASTRAN release, an attempt is made to avoid the above dilemma. The 
orthogonality convergence criterion is based on a ratio instead of the finite difference of two very 
close numbers. However, presently there is not enough test data to verify that is a good fix. 


NEW LOGIC FOR MATRIX TRANSPOSE 

Matrix transpose of a matrix which is too big to reside completely in the computer core 
memory is not an easy task. It can be done, but it may use up lots of cpu time. Lots of I/O may 
be involved here, and perhaps a high percentage of system paging if a virtual machine is used. 
The problem here is how to do the matrix transpose in seconds instead of minutes, or in minutes 
instead of hours. A good algorithm here is a treasure, and quite often, it becomes a "company 
proprietary" product. 

The out-of-core matrix transpose in NASTRAN is handled by the subroutine TRNSP. The 
algorithm there is amazingly powerful. The only drawback is that it uses up to nine scratch files. 
The scratch files are supplied by the calling routines. The more scratch files passing over to 
TRNSP, the bigger matrix transpose TRNSP can do. There is no check in TRNSP of the actual 
scratch files requirement. Again, there are lots of I/O, data packing, and unpacking involved. 

If only l/50th of the matrix can be loaded into the computer core memory space at one 
time, TRNSP will complete the transpose task in 50 passes. A new matrix transpose subroutine 
TRNSPS has been written for 92 COSMIC/NASTRAN. If the passes exceed seven, TRNSP will 
switch over to TRNSPS automatically, if and only if DIAG 41 is not turned on by the user. 
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TRNSPS uses only one scratch file. The I/O department and the data packing and unpacking are 
greatly reduced. The new TRNSPS is two to four times faster than the original TRNSP. 


MPYAD, MPY4T, AND MPYDRI 

Matrix multiplication and addition are basically the most important and most widely used 
tools for a finite element program. If the matrices are small and can reside completely in core, 
multiply-add has no problem; and three simple do-loops will complete the job. Again, if the 
matrices are bigger than the computer core can hold, matrix multiplication and addition have to 
be done by parts, and there will be many passes, many I/O operations, and many row- and 
column-packings and unpackings. The problem can become I/O bound, and lots of cpu time will 
be spent on getting and saving the intermediate results. The situation is further complicated in 
NASTRAN in that the matrices can be of different types (single precision or double precision; 
real or complex), or of different forms (rectangular, square, diagonal, identity that may or may 
not exit, lower or upper triangular, row vector, or symmetric), or transpose and non-transpose 
matrices. 

There are five multiply-add (MPYAD) methods in NASTRAN, two methods for the 
non-transpose case (MPY1NT and MPY2NT) and three for MPYAD with transpose (MPY1T, 
MPY2T and MPY3T). NASTRAN selects internally the best method to use based on the cpu 
time requirement of each method that fits best for a given matrices-and-core environment TTie 
cpu time requirement is a function of size and shape (rows and columns), form and density of 
the matrices, and core space. The cpu time requirement is also a function of the relative sizes and 
shapes of the matrices. That is, matrix A may be very big and cannot reside in core, while matrix 
B is small, and matrix B may or may not be loaded entirely into the available core space. Or vice 
versa. All five MPYAD methods are written in assembly languages for four out of five COSMIC 
supported machines, except VAX. 

Matrix transpose is supposedly very slow. The programmer manual recommends matrix 
transpose be done via MPYAD with transpose, and an identity matrix. This turns out not to be 
very efficient. Two test matrices A and B, 5166x5166 each, using the best of MPY1NT or 
MPY2NT, can be multiplied together in 470 cpu seconds (on the VAX machine), while 
A-transpose times B, using the best of MPY1T, MPY2T, and MPY3T requires 6680 cpu seconds. 
That is 12 times longer! If matrix A is transposed first by TRNSP routine (TRNSPS is not used), 
then followed by MPYAD without transposing, the total cpu time can be cut to half. The only 
problem here is that at least one extra scratch file is needed to save the transpose file, and in 
many instances, there is no extra scratch file available. 

MPY1NT and MPY1T share much common logic, and operate quite similarly. The same 
holds true for the MPY2NT and MPY2T pair. MPY3T is a third method for the transpose case. 
The best method for the test matrices A and B above, with transpose, is MPY2T. One would 
think that since NASTRAN stores a matrix by column, and that a column of matrix A transpose 
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is a row of A, the row-and-column multiplication (matrix A in row and matrix B in column) 
should be very fast, very smooth, and very convenient. One could almost feel and touch the 
natural flow of the multiplication algorithm. But what one can feel or imagine, is not what a 
computer sees. The row-and-column multiplication produces only one element in the resulting 
matrix. There are more than 26 million (5166 2 ) double precision elements to go. Anyway, the fact 
is that MPY2T takes 12 times longer than MPY2NT. 

Several methods, several logics, and several new algorithms were developed and tested 
to beat the clock set by MPY2T. The ultimate goal is to match the MPY2NT performance if 
possible. Finally, after many trials, a fourth method, MPY4T, was developed based on the scheme 
similar to MPY2T, except that matrix B (and matrix C, to be added if it is present), and the 
resulting matrix D are processed by column instead of by element. Of course, the logic 
MPY4T becomes much more complicated and the open core must be rearranged. But MPY4T 
is three to five times faster than MPY2T. In the 92 NASTRAN release, MPY4T will be 
automatically substituted for MPY2T, unless DIAG 41 is turned on by the user. MPY4T is 
written in FORTRAN, and it is machine independent. 

The original MPYAD does not take advantage of certain types of matrices. For example, 
the transpose of matrix A which is symmetric, need not go through the transpose route. (This is 
already checked in the 91 release.) A new subroutine, MPYDRI, is added to the 92 release to 
handle special cases involving diagonal matrix, row vector, and identity matrix. 

The correct handling of these special matrices expedites the MPYAD process by 
manyfold. 


NEXT IMPROVEMENTS UNDER CONSIDERATION 

The internal improvements in the 1992 NASTRAN release tackle a few important, and 
often-used, basic finite element tools with satisfactory results. However, there are many more 
areas in NASTRAN that can be explored. There are still several areas involving FBS that have 
not been touched. The matrix decomposition process could be improved. All the complex 
computations involving complex FBS, complex decomposition, complex FEER method, and 
more, are targets for the next improvements. Nevertheless, the internal improvements in 1992 
NASTRAN release represent the beginning of an extraordinary effort to bring NASTRAN up to 
par. 
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APPENDIX A 


Demo problem D03012A was used in most of the following tests. The D03012A demo 
produced a double precision KGGX matrix of size (5166 x 5166) and a KAA (2380 x 2380), 
double precision. The trailer of the KGGX matrix in some cases had to changed from 
"symmetric" to "square" so that NASTRAN did not take the symmetric route processing the 
modules under investigation. Tests were done on a WAX/VMS machine, unless stated otherwise. 
In most cases, HICORE is 350,000 words. 


FBS test: KAA 


1991 Version 

1992 First Version 

1992 Second Version 

5644 cpu seconds 

3508 cpu seconds 

3120 cpu seconds 


FEER method 


1991 VAX Version 

1992 VAX Version 

GINO Improvement 

Open Core Not Used 

Open Core Used 

1043 cpu seconds 

978 cpu seconds 

907 cpu seconds 

763 cpu seconds 

1991 CRAY Version 

1991 CRAY Version Plus Changes - Open Core Used 

45.7 cpu seconds 

21.8 cpu seconds 


MPYAD, KGGX(transpose) * KGGX + KGGX 


DIAG 41 On 

With MYADT 
(Obsolete) 

MPYAD With 
New TRNSPS 

MPYAD With 
New MPY4T 

If Symmetric 
Matrix Allowed 

6681 cpu secs. 

3871 cpu secs. 

21 14 cpu secs. 

1358 cpu secs. 

469 cpu secs. 
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EVOLUTION OF A NASTRAN TRAINER 

H.R. Grooms, P.J. Hinz, and M.A. Collier 

Rockwell International 
Downey, California 


INTRODUCTION 

This paper traces the development of a NASTRAN training system. It encompasses the design and 
organization of the program, including the static and dynamic modules. A discussion of how user 
feedback, in the form of questionnaire responses, was used to evaluate and improve the trainer is included. 


BACKGROUND 


The user-friendly NASTRAN trainer was originally designed as a segment of a larger system 
(ref. 1). After the static module was developed, used, and evaluated (ref. 2), it became clear that the trainer 
concept readily lent itself to a wide range of applications. 


The NASTRAN trainer (figure 1) was initially conceived as a low-cost, convenient tool for giving 
engineers who were novices in finite element analysis a few practical applications of the method 
Although several very good short courses and classes on NASTRAN are available, most are offered 
periodically and cost $200 to $800 per person. These classes usually require the engineer to set aside his 
current work assignment and devote his full time to the class for anywhere from a couple of days to a 
couple of weeks. When funds are low or schedules are tight, the money or time required for these classes 

can be insurmountable barriers. 


Various researchers have developed computer programs for structural analysis and design 
applications. Ginsburg (ref. 3) addresses computer literacy, and Woodward and Moms dis cuss improved 
productivity through interactive processing (ref. 4). Wilson and Holt (ref. 5) developed a system o 
computer-assisted learning in structural engineering. Sadd and Rolph (ref. 6) describe the various ways in 
which design engineers could be trained to use the finite element method. Self-adapting menus for 
computer-aided design (CAD) software are covered by Ginsburg (ref. 7). 

Bykat (ref. 8) is developing a system that will have features for training, analysis control, and 
interrogation. 


STATIC MODULE 


The static module was developed to provide the user a variety of different types of problems. The ten 
problems generally increase in difficulty as their numbers increase. Table I describes each problem, an 
figures 2 and 3 show the problems. Figure 4 gives the classical solution for static example 6. 


The user may work the problems in any order. 
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DYNAMIC MODULE 


The dynamic module contains eight problems (table II) of increasing complexity. The problems are 
shown in figures 5 and 6. A classical solution for dynamic problem 4 is presented in figure 7. The 
examples were selected to allow the user to decide on: 

1 . Grid fineness 

2. Mass representation 

3. Number of degrees of freedom retained 

4. Particular degrees of freedom retained 

Each of these decisions can have a significant bearing on the accuracy of the eigensolution. 

A complete description of this module is given in ref. 9. 


TRAINER ORGANIZATION 

The NASTRAN Environment (NE) is written for the IBM computer system running MVS/ESA 
SP 4.1.0 using TSO/E 2.1.0. It uses the features provided by the dialogue management services under 
ISPF/PDF to the fullest extent. This includes panels, skeletons, CLISTs, and tutorial services. In addition, 
the trainer requires VS/FORTRAN and VS/Pascal compilers if the executable code is not directly 
portable. Newer versions of any of these services should not invalidate the NE if the improvements are 
upwardly compatible. 

The NE has job setups to execute MSC/NASTRAN and COSMIC/NASTRAN on IBM and/or 
MSC/NASTRAN on the Cray running Unicos 6.1. At least one of these programs must be available for 
the NE to be used as intended. These codes are not delivered with the NE. 

The NE comprises 15 datasets; 14 are partitioned and 1 is sequential. These datasets are listed below 
with a brief explanation of their contents. Figure 8 illustrates the organization of the NE datasets. 

• ALTER— Rigid format alter library for NASTRAN. Must be updated with each release of a new 
version of NASTRAN. Not a requirement for the trainer and most NASTRAN users. 

• CLIST — Procedural commands to invoke the NE, allocate and manage datasets, create and submit 
batch jobs, invoke the SPF editor and initialize profile variables. 

• MSGS — Messages that appear on panels for information, caution, and warning. 

• PNLZ (Environment) — Panels that serve as the user interface. All information that the user inputs 
and the system outputs is through panels. This includes data entry screens, system news, help 
sections, user manuals, problem descriptions, and classical solutions. 

• DOC — Documentation on efforts to develop NASTRAN expert systems. 

• FORT — FORTRAN programs and subroutines that compute the classical solutions and extract 
system jobcard information. 
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. INP Input decks that are example solutions to the NASTRAN trainer problems. 

• LOAD — Load modules of the compiled and linked FORTRAN and Pascal programs. 

. LOG Log of NASTRAN trainer usage for each NASTRAN trainer job submitted. 

• OBJ — Object modules of the compiled FORTRAN and Pascal programs. 

. OUT Output decks that are example solutions to the NASTRAN trainer problems. 

. pvs — Pascal program that produces a report of the usage of the NASTRAN trainer. 


FUTURE ENHANCEMENTS AND USER FEEDBACK 


Modules for elastic stability (buckling) and substructuring are in the Panning! stage s jJ dltlons 

are planned as self-contained units that can be used by anyone who has completed the static modu . 


Users of the static module were asked to fill out a questionnaire. The questions 
shown in figure 9. Another questionnaire is currently being used to solicit opinions 

module. 


and responses are 
about the dynamic 


CONCLUSIONS 


The NASTRAN trainer has been used by a number of engineers, who found it to be a versatile low- 
cost tool It is particularly helpful in bridging the gap from theory to practical application of the finite 
elemcm method for straccural analysis. The program, along with documen.at.on, ,s ava.lable through 


COSMIC. 
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TABLE I. -STATIC EXAMPLE PROBLEMS 


Example 

Description 

Significant Features 

1 

Statically determinate plane truss subjected to point load 

Bar elements, stability constraints 

2 

Beam simply supported on one end and fixed at the other 
subjected to point load 

Beam elements 

3 

Beam fixed at both ends subjected to through-the-depth 
temperature difference 

Temperature input 

4 

Plane frame subjected to point load 

Half-model, symmetric, and antisymmetric loads 

5 

Simply supported beam subjected to temperature pattern 

Half-model, temperature distribution decomposed into 
symmetric and antisymmetric parts 

6 

Plate with hole in center subjected to in-plane load 

Plane stress, quarter- model, fine grid around hole 

7 

Simply supported square plate subjected to out-of-plane 
point load at center 

Plate-bending elements, quarter-model 

8 

Three-dimensional frame subjected to point load 

Tapered beams, three-dimensional 

9 

Cylindrical shell subjected to hydrostatic loading 

Three-dimensional simulation of curved surface using 
flat elements 

10 

Cylindrical shell with ring frames closed at both ends 
subjected to internal pressure 

Self-equilibrating loading, three-dimensional 


TABLE II. -DYNAMIC EXAMPLE PROBLEMS 


Example 

Description 

Significant Features 

1 

Beam simply supported on both ends with lumped mass in 
middle 

Motion in one plane only, lumped mass only 

2 

Beam simply supported on both ends with uniformly 
distributed mass 

Motion in one plane only, distributed mass 

3 

Beam fixed on one end with a lumped mass at the free end 

Motion in any direction, lumped mass only 

4 

Beam fixed on one end with a uniformly distributed mass 

Motion in any direction with uniformly distributed mass 

5 

Rectangular plate clamped on one edge, all other edges free, 
with a uniformly distributed mass 

Plate bending with distributed mass 

6 

Rectangular plate, free-free with uniformly distributed mass 

Free-free (implies six modes with zero frequency) 

7 

Two beams connected by springs, each with distributed and 
lumped mass 

Multibody problem, free-free 

8 

Problem 7 with a forcing function added 

Forcing function 
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Static Problem 3 - Beam Fixed at Both Ends 


Static Problem 1 - Two-Dimensional Truss with Temperature Loading 



FIGURE 2.-STATIC PROBLEMS 1 TO 4 
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Static Problem 6 - Plate With Hole in Center 



o-j-= Stress at Distance From Hole (psi) ^ pp 

= Stress (X) at Hole (psi) Where o nom - t p . 2r ) 

o Y = Stress (Y) at Hole (psi) k = 3 00 - 3.13(-^) + 3.66(-|-) 2 - 1.53(^) 3 

Calculate ay at Point b: 

P 

oy = Stress at Distance From Hole = -y- 

Oy =0 b =-Oy 

Reference: Roark and Young. Formulas tor Stress and Strain, 5th Edition, p. 594. ^td 920210-3135 


FIGURE 4 .-CLASSICAL SOLUTION FOR STATIC PROBLEM 6 
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FIGURE 5.-DYNAMIC PROBLEMS 1 TO 4 
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MTD 920211-3137 


FIGURE 6-DYNAMIC PROBLEMS 5 TO 8 
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Dynamic Problem 4 - Cantilever Beam With 
Uniformly Distributed Mass 



w= ib/in 
A, E 



Cross Section 


E = Modulus of Elasticity (psi) 

I = Moment of Inertia (in. 4 ) 

p = Distributed Mass (lb-sec 2 /in. 2 ) 
L = Length of Beam (in.) 


<o n = Natural Frequency - Angular (rad/sec) 

f n = Natural Frequency (cycles/sec or hertz) 

Reference: Flugge, W.: Handbook of Engineering 
Mechanics. McGraw-Hill 
1962, pp. 61-8. 


Classical Solution: 

Calculate Natural Frequencies: 




f 



(Fundamental Mode) 

(Higher Order Modes) 
(n > 1) 


Where n =1,2,... (the Mode Number) 


MTO 920212-3138 


FIGURE 7.-CLASSICAL SOLUTION FOR DYNAMIC PROBLEM 4 
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FIGURE 8.-NASTRAN ENVIRONMENT DATASET ORGANIZATION 
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Critique of NASTRAN Trainer 


1 Was using this system a worthwhile expenditure of your time? 

a. Yes (89%) 

b. No (0%) 

c. Undecided (11%) 

2 How much total time would you estimate that you spent using the trainer? 

60 hours 

3 How much total time would you have spent (estimate) to gain this knowledge if the trainer had not been 
available 135 hours 

4 The number of examples was 

a. Too few (17%) 

b. Too many (6%) 

c. About right (77%) 

5 The system was 

a. Too simple (17%) 

b. Too complicated (6%) 

c. About right (77%) 

6 Could the trainer be improved by adding other topics? 

a. Yes (67%) 

b. No (22%) 

c. Maybe (11%) 

7 Which section, if any, should be expanded upon? 

(Wide variety of responses.) 

8 How often (average) did you invoke the NASTRAN documentation manual section? 

a. Never (44%) 

b. 0-2 times/example (22%) 

c. More than 2 times/example (34%) 

9 Was the NASTRAN documentation section useful? 

a. Yes (38%) 

b. No (33%) 

c. Never used it (29%) 

10 How often did you use (average) the printed COSMIC or MSC NASTRAN manuals? 

a. Never (6%) 

b. 0-2 times/example (17%) 

c. More than 2 times/example (77%) 

1 1 Please add any additional comments you desire. 

(Responses vary from "great" to "give us more advanced problems.") 
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ANIMATION OF FINITE ELEMENT MODELS AND RESULTS 


Robert R. Lipman 


David Taylor Research Center 

Computational Signatures and Structures Branch (Code 1282) 
Bethesda, Maryland 20084-5000 


SUMMARY 

Several years ago, the phrase ‘visualization in scientific computing’ (ref. 1) was coined for 
what we used to call computer graphics. Although computer graphics is part of visualization, 
visualization encompasses computer graphics hardware and software, network communications, 
user interfaces, computer-aided-design, and more. The purpose of visualization is to provide 
insight into the engineer’s models and calculations. Animation of finite element models and 
results is a visualization process that can provide the insight. 

The paper is not intended to be a complete review of computer hardware and software 
that can be used for animation of finite element model and results, but is instead a 
demonstration of the benefits of visualization using selected hardware and software. Opinions 
expressed are solely those of the author and are not those of the David Taylor Research Center, 
the Navy, or the Department of Defense. Good reviews of visualization hardware and software 
can be found in the following journals: Computer Graphics World, Supercomputing Review, 
IEEE Computer Graphics and Applications, and CAE Computer-Aided Engineering. A 
videotape showing visualization and animation of finite element models and results is an integral 
part of this paper although it is not included in the proceedings. 


INTRODUCTION 


Visualization and animation give an engineer insight into his finite element model and 
results. Wire-frame plots of a finite element mesh do not convey the sense that a ‘real’ structure 
has been modeled. We do not live in a wire-frame world. We live in a world of color, light, 
shading, and perspective. A beam is not a line between two points. A beam has a web and a 
flange of substantial size and cross-section. The transient motion of a structure cannot be 
determined from static plots at selected time steps or plots of the response of a node versus 
time. Visualization and animation can be used to show an engineer the realistic configuration 
and response of a structure. 

The earliest animations of finite element analysis results were made by painstakingly 
recording a sequence of static plots on film. Some of the first computer animations of finite 
element analysis results were made on Evans & Sutherland graphics hardware (ref. 2). Today, 
with the price/performance of computer graphics hardware so low and visualization software 
packages becoming more mature, finite element model visualization and animation is now a 
desktop tool. 
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HARDWARE 


It seems that computer workstation vendors are announcing faster and less expensive 
hardware almost every day. Current Cray X-MP supercomputer computational speeds should be 
available in desktop computer workstations in 1-2 years. By the time this paper is presented, 
the 1-2 year time frame may have been reduced from years to months. The cost of memory and 
hard disk storage is also falling. 

As important as raw computational power is to animation and visualization, graphics 
speed is equally important. Graphics speed, usually quoted in polygons per second, is not 
increasing as fast as computational speed. Rather, the cost of current graphics power is getting 
less expensive. Current peak graphics speeds of 200,000 polygons per second can be found on 
the top-of-the-line graphics workstations. The user must beware of the type of polygon that the 
vendor uses when quoting graphics speed. Quotes of one million polygons per second arc usually 
for highly optimized meshes of triangles without light-source shading. For animation of finite ele- 
ment models, the graphics speed for independent quadrilateral polygons is more relevant. 

When computation and visualization take place in a distributed environment, communica- 
tions speed between the client and server is another important issue. Today, animations of finite 
clement analysis results are typically done in a batch mode. The analysis is done on a large 
mainframe or supercomputer and the results are sent to a workstation to be used with visualiza- 
tion software. In the future, the two processes of analysis and visualization will be more tightly 
coupled where the analysis and visualization are being computed concurrently. For this scenario 
to take place, much higher network communications speed between the computational server 
and visualization server will be necessary than current local- and wide-area networks provide. 

Finally, animation sequences have to be recorded to videotape. There are two methods 
for recording computer graphics animations on videotape: real-time and frame-by-frame. For 
real-time recording, the computer graphics display is converted, in real-time, to a television sig- 
nal suitable for recording on videotape and being displayed on a regular television monitor. 
Therefore, whatever is being displayed on the computer graphics display can be recorded to 
videotape. If the graphics speed is fast enough to animate a finite element model in real-time, 
then this process is sufficient. 

With graphics hardware that is not fast enough and with visualization software that has a 
rendering capability, then frame-by-frame recording can be used. The visualization software 
renders individual images that are recorded one-by-one on videotape. The result is a continuous 
animation sequence. Frame-by-frame recording also produces higher quality animations and 
renderings than real-time animation. 


SOFTWARE 

Visualization software packages can be separated into three categories: general-purpose, 
modular, and application specific. The types of data that the packages can visualize are usually 
either structured or unstructured grids. A structured grid is typical for finite-difference applica- 
tions such as computational fluid dynamics. A finite element mesh is an example of an unstruc- 
tured grid. Many of the general-purpose visualization packages (PV-Wave, Spyglass, Data 
Visualizer) are very good with structured grids and less useful for finite element applications. 
There are also several application specific visualization packages (Fieldview, Fast, Plot3D) that 
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can be used with only structured grids. The modular visualization packages (A VS, Iris Explorer, 
apE) allow users to write their own applications specific software modules to be integrated with 
the visualization software. 

There are few choices for application-specific visualization packages for finite element 
analysis animation. The popular finite element pre and postprocessors (Patran, I-DEAS) have 
animation capabilities, but are not oriented to the visualization process and to recording video- 
tapes. FOTO (ref. 3) is a data visualization software package geared towards finite clement 
models and animation. FOTO was used to make the videotape that is part of this paper. FOTO 
is easy to interface with analysis codes; is user-definable menu driven; has many visualization 
types including: color, displacement, contour lines, vectors, transparency, and culling operators; 
and has a tightly coupled videotape system. 

There are also free visualization software packages available from the national supercom- 
puter centers. Because they are free, the source code is provided allowing the user to tailor the 
code to his application. However, because they are free, the user will not get the same type of 
support or updates to the software that a commercially available package would provide. 


THE FUTURE 

What is currently possible for finite element animation and visualization is not the final 
product, but only a step towards a more interactive, dynamic environment for doing analysis and 
visualization. In the future, analysis and visualization will occur concurrently in near real-time 
and the engineer will have the capability to interact with the analysis by changing the finite ele- 
ment model as the computations are taking place to explore new configurations ol the model. 
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ACCURACY OF THE TRIA3 THICK SHELL ELEMENT 

William. R. Case, Marco Concha and Mark McGinnis 
NASA/Goddard Space Flight Center 


SUMMARY 


The accuracy of the new TRIA3 thick shell element is assessed via comparison with 
a theoretical solution for thick homogeneous and honeycomb flat simply supported 
plates under the action of a uniform pressure load. The theoretical thick pLate 
solution is based on the theory developed by Reissner and includes the effects of 
transverse shear flexibility which are not included in the thin plate solutions based 
on Kirchoff plate theory. In addition, the TRIA3 is assessed using a set of finite 
element test problems developed by the MacNeal-Schwendler Corp (MSC) 
Comparison of the COSMIC TRIA3 element as well as those from MSC and 
Universal Analytics Inc. (UAI), for these test problems is presented. The current 
COSMIC TRIA3 element is shown to have excellent comparison with both the 
theo reti cal solutions and also those from the two commercial versions of 
NASTRAN with which it was compared. 


INTRODUCTION 


The TRIA3 thick shell element was added to the 1990 release of COSMIC 
NASTRAN. Along with the QUAD4, the two new shell elements represent a 
significant increase in the capability of COSMIC NASTRAN to model complicated 
shell structures. The deficiencies of the original TRIA1.2 and QUAD1.2 shell 
elements have been recognized for years and have been reported in the literature. 
At the Goddard Space Flight Center (GSFC), the triangular and quadrilateral shell 
elements are used in virtually all structural analyses of our spacecraft and related 
hardware. Typical applications are for the modeling of cylindrical shells and flat 
plates made of honeycomb or machined, lightweighted, metal that make up the 
structure of spacecraft and scientific instruments. In some cases these models 
require that the effects of transverse shear flexibility be included due to their 
thickness. The TRIA3 and QUAD4 elements include these effects. The QUAD4 
element has, in addition, an improved membrane capability for in-plane loading. 
The TRIA3 element, due to it's limited number of degrees of freedom retains the 
constant strain membrane capability of the older TRIA1 and TRIA2 elements. This 
necessitates finer meshes for in plane loading cases than would be required when 
using the QUAD4 element. 

The purpose of the study reported herein is to assess the accuracy of the TRIA3 
element in modeling a variety of situations involving both solid cross-section plates 
as well as those constructed of honeycomb. An identical study for the QUAD4 
element was reported in the 18th NASTRAN User's Colloquium and is documented 
m reference 1. As with the QUAD4 study, the three goals of the TRIA3 study were 
to determine: 
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a) what is the rate of convergence to the theoretical solution as the mesh is 
refined 

b) whether the element exhibits sensitivity to aspect ratios significantly 
different than 1.0 

c) how the element behaves in a wide variety of modeling situations, such as 
those included in the MSC element test library (discussed below). 

The first two questions were addressed in the same manner as several other studies 
reported by one of the authors in prior NASTRAN colloquia (references 1 - 3). The 
procedure used in those studies, and followed here also, is to isolate the effects of 
mesh refinement and aspect ratio. That is, the mesh refinement study is done using 
elements with an aspect ratio of 1.0. Then, once a fine enough mesh has been 
reached such that the errors are small, the effects of aspect ratio can be investigated 
by keeping the mesh the same (i.e. same number of elements) and varying the 
overall dimensions of the problem, thus resulting in each element aspect ratio 
changing. Obviously, in order to accomplish this latter step there must be a 
theoretical solution (or some other equally acceptable comparison solution) to the 
problem with which to compare the finite element model results. This is needed 
since, at each step, a problem of different dimensions (and therefore different 
theoretical solution) is being modeled. 

The above tests are important in that they show the rate of convergence toward the 
theoretical solution as the mesh is refined. Those tests, however, are not sufficient 
to completely test the accuracy of a finite element since they do not test irregular 
geometries, or a variety of loadings or material properties. The MSC has developed 
a comprehensive set of problems for testing finite elements in a variety of situations 
(reference 4). The library of problems consists of 15 test problems for shell 
elements that cover all of the parameters mentioned above. This element test 
library was used to test the TRIA3 element as was done for the QUAD4 element 
reported in reference 1. 


RESULTS OF MESH AND ASPECT RATIO STUDY 


For the mesh and aspect ratio study a theoretical comparison solution is hig hly 
desirable. Since the effects of transverse shear flexibility are included in the TRIA3 
element formulation, a theoretical solution for moderately thick plates, based on 
Reissner (or Mindlin) thick plate theory is also desirable. Such a solution is given in 
references 5 and 6 for rectangular simply supported thick plates under the action of 
a pressure load. Thus, this problem was used for the mesh and aspect ratio portions 
ofthe study. 

Figure 1 defines the geometry, coordinate system, boundary conditions and loading 
for the rectangular plate. The thickness indicates a moderately thick plate of length 
to thickness ratio of 20. The effect of transverse shear flexibility is only 
approximately 1% on the maximum displacement but is important in discerning the 
quality of the convergence of the finite element results to tne exact theoretical 
solution. By exact is meant the theoretical basis for the TRIA3 element, which is 
expressed in the Reissner thick plate theory. Figure 2 shows the finite element 
mesh geometry used in the mesh and aspect ratio studies. Due to symmetry only 
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one quarter of the plate was modeled. The 4x4 mesh shown in figure 2 is an 
example only; the mesh was varied during the mesh study. However, as was done 
for all problems, the quad areas were subdivided into triangles in the alternating 
orientation shown in figure 2. 

Figures 3a - 3c show characteristics of the theoretical solution. As indicated in 
figure 3a the central displacement solution is represented as an infinite series of 
hyperbolic functions. A FORTRAN computer program was written to compute the 
theoretical solutions for displacements (using the series shown) as well as stresses 
(solution not shown). Figures 3b and 3c show the stiffness parameters needed in the 
theoretical solution for the homogeneous and honeycomb plates. For the 
honeycomb plate, two different core stiffnesses were investigated. The stiffer one is 
representative of aluminum honeycomb construction that has been used at the 
GSFC. The more flexible one was chosen because it represents a core flexibility 
that is quite low and was expected to be a more critical test of the TRIA3’s shear 
flexibility formulation. 

The results of the mesh study, showing the convergence of the TRIA3 solutions to 
the theoretical, are presented in tabular form in tables 1-2 and in graphical form in 
figures 4 - 7. Both formats show percent error in displacement at the center of the 
plate as a function of mesh refinement. Results are included for COSMIC 9.0, UAI 
11.1 and MSC 66A NASTRAN. The tables merely give exact numbers (along with 
the theoretical displacements) and the figures contain the same error information, 
but in graphic form. Figures 4 and 5 and table 1 are the results for the 
homogeneous plate. The difference between the results in figures 4 and 5 (and that 
in the two parts of table 1) is that figure 4 (and the top half of table 1) is for a 
solution in which shear flexibility is included and figure 5 (and the bottom half of 
table 1) neglects shear flexibility. These two situations were investigated to test the 
MID3 option on the PSHELL NASTRAN bulk data deck card which allows the 
effects of shear flexibility to be ignored if MID3 is left blank. As seen in figures 4 
and 5 the NASTRAN results converge very rapidly with mesh refinement tor 
COSMIC 9.0, MSC 66 A and UAI 11 . 1 . As seen, all versions converge to less than 
1% error for a mesh size of 8 x 8. 

Figures 6 and 7 and table 2 are the results for the honeycomb plate. Figure 6 (and 
the top half of table 2) are for the honeycomb plate with the stiffer core and figure 7 
(and the botto m h alf of table 2) are for the more flexible core. As seen in figures 6 
and 7 the NASTRAN results for COSMIC 9.0 and the two commercial NASTRAN 
versions converge very rapidly for the two honeycomb plates as they did for the 
homogeneous plate. 

In order to test the TRIA3's sensitivity to aspect ratio, the model with a 12 x 12 
mesh was run in which the plate side dimension in the x direction was varied. This 
causes the element aspect ratio to vary while maintaining a constant mesh in an 
attempt to prevent mesh refinement errors from significantly affecting the results. 

As seen in tables 1 and 2, the TRIA3 results with a 12 x 12 mesh (and aspect ratio of 
1.0) have very little error. The results of the aspect ratio study are presented in 
figures 8-10 and tables 3-5. Tables 3-5 give percent error in the displacement at 
the center of the plate versus aspect ratio for a model with a mesh of 12 x 12 TRIA3 
elements (over one quarter of tne plate). As mentioned above, the aspect ratio was 
varied by changing the dimension of the plate along the x axis. For example, the 
results for the aspect ratio of 10 are for a plate (and all TRIA3 elements) that is 10 
times as long in the x direction as in the y direction. Therefore, the theoretical 
solution changes with aspect ratio. Figure 8 and table 3 are for the homogeneous 
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plate (with transverse shear flexibility) while figure 9 and table 4 are for the stiff 
core honeycomb plate and figure 10 and table 5 are for the more flexible core 
honeycomb plate. Investigation of the percent error in the tables, as well. as in 
figures 8 - 10 show that the TRIA3 has essentially no aspect ratio sensitivity over the 
range investigated. 

Based on the above results, the COSMIC TRIA3 element is seen to give very 
accurate results for the displacements in the problem investigated, both in 
comparison to the exact theory and in comparison to the two commercial versions of 
NASTRAN that we have at the GSFC. Although the results are not presented 
herein, similarly accurate results were obtained for the shear and moment stress 
resultants as well. In addition, the rates of convergence for the TRIA3 compare 
quite favorably with that found for the QUAD4 in reference 1 for this plate bending 
problem. 


RESULTS OF TESTING USING THE MSC ELEMENT TEST LIBRARY 


As mentioned earlier, the mesh and aspect ratio studies, while a verv useful tool in 
the evaluation of an element, do not test all of the important variables that affect 
accuracy in a finite element solution. The MSC element test library mentioned 
above represents a rather exhaustive series of tests that include many of the element 
related parameters which affect the accuracy of a finite element solution. 

Reference 4 gives a detailed description of the test problems along with theoretical 
answers and the results of the testing on several MSC elements. The reader should 
consult reference 4 for a complete description of the various problems in the test 
series. The portion of this series of element tests that relate to shell elements was 
run by the authors on the TRIA3 elements contained in COSMIC 9.0, UAI 11.1 and 
MSC 66A. As the MSC does in their report, the results are presented in detail and 
also in a summary form in which the element is given a letter grade of A through F 
based on the magnitude of the error. Table 6 shows the summary results for the 15 
tests in the series ranging from a simple patch test to modeling of beams (using the 
TRIA3 element through the depth) and various plates and shells. The meaning of 
the letter grades is given at the bottom of the table. As pointed out in reference 4, a 
failing grade for an element in one test is not a reason to dismiss the element. For 
one thing, the test scores would improve with mesh refinement; the mesh used in 
most of the problems was quite coarse. Of importance in this discussion is not the 
actual grades listed in table 6 but the comparison of the COSMIC grades with those 
from the other two programs. As seen in table 6, the COSMIC TRIA3 element is as 
good as, or better than, those of the commercial programs. All of the low marks (D 
or F) are apparently due to the constant strain membrane portion of the TRIA3 
element and the low order mesh used in those problems. For example, the straight 
beam bending, with in-plane loading, had only one TRIA3 through the thickness. 
This was done to keep the same mesh as MSC used for the QUAD4 element tests, 
and was also done in reference 1. Refining the mesh would have improved the 
answers to any degree of accuracy desired; the low grades are not indicative of any 
failure of the element to converge. Although not shown in table 6, the old TRIA2 
element (included in reference 4) has a D or F grade in 9 of the 15 problems. The 
twisted beam test (number 11 in table 6) is really used to test the effect of warp on 
quadrilateral elements, which is not applicable for the TRIA3 element. 
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CONCLUSIONS 


The COSMIC TRIA3 general purpose flat shell element has been shown to be an 
excellent element and, together with the QUAD4 quadrilateral flat shell element 
significantly enhances the usefulness of COSMIC NASTRAN. The element has ’ 
been shown to compare excellently with those available in two commercial versions 
of NASTRAN that are currently being used at the GSFC. 
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List of Symbols 


w = plate displacement in the z direction 
x,y,z = coordinate directions 
p = pressure load on the plate in the z direction 
a, b = plate dimensions (length, width) 
t = overall plate thickness 
D = plate bending rigidity (see Figures 3b, 3c) 

Cs, C n = plate shear stiffness (see Figures 3b, 3c) 

tf = thickness of face sheets for honeycomb plate 

tc = thickness of the core for honeycomb plate 

N x = number of elements in x direction in one quarter of plate 

Ny = number of elements in y direction in one quarter of plate 

AR e = element aspect ratio (see Figure 2) 

E = Young's modulus 
G = shear modulus 
v = Poisson's Ratio 
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TABLE 1: TRIA3 Error in Displacement at Center of Plate 

Mesh Size Study (Element Aspect Ratio 1.0) 
Simply-Supported, Homogeneous Plate Under Uniform Pressure Load 


Theoretical Displacements 
With Transverse Shear Flexibility: 3.571 x 10' 5 m 

(1.406 x 10'3 in.) 

Without Transverse Shear Flexibility: 3.529 x 10*5 m 

(1.390 x 10-3 in.) 


% Error 

Cosmic UAI MSC 

Mesh 90 Ver. 11.1 Ver. 66A 

With Transverse Shear Flexibility 


lxl 

39.33 

27.64 

16.62 

2x2 

13.63 

11.36 

9.01 

4x4 

3.29 

2.77 

2.06 

8x8 

0.01 

0.55 

0.34 

12x12 

0.00 

0.13 

0.04 

Without Transverse Shear Flexibility 
lxl 40.56 28.37 

17.45 

2x2 

14.31 

11.72 

9.52 

4x4 

3.74 

3.01 

2.43 

8x8 

0.95 

0.76 

0.62 

12x12 

0.42 

0.34 

0.27 
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TABLE 2: TRIA3 Error in Displacement at Center of Plate 

Mesh Size Study (Element Aspect Ratio 1.0) 
Simply-Supported, Honeycomb Plate 
Under Uniform Pressure Load with Transverse Shear Flexibility 

Theoretical Displacements 
G z = 1 .5 1 7x 1 08 N/m2 : 2.422x 1 0‘3 m 

(9.535x10"2 in.) 

G z = 1.379x107 N/m2 ; 3.102x10-3 m 

(1.221x10-1 in.) 


Mesh 

Cosmic 

90 

% Error 
UAI 

Ver. 11.1 

MSC 
Ver. 66A 

Gz = 

1.517x108 N/m2 (22000 psi) 


lxl 

38.28 

27.13 

16.08 

2x2 

13.36 

11.36 

8.88 

4x4 

3.35 

2.92 

2.16 

8x8 

0.81 

0.73 

0.52 

12x12 

: 0.37 

0.33 

0.24 

G z = 

1.379x107 N/m2 (2000 psi) 


lxl 

24.07 

17.82 

7.37 

2x2 

9.71 

8.83 

6.35 

4x4 

2.48 

2.26 

1.60 

8x8 

0.60 

0.55 

0.38 

12x12 

0.30 

0.28 

0.02 
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TABLE 3: TRIA3 Error in Displacement at Center of Plate 

Aspect Ratio Study (12 x 12 Mesh) 

Homogeneous, Simply-Supported Plate 
Under Uniform Pressure Load with Transverse Shear Flexibility 


AR 

theoretical w, 
m (in.) 

Cosmic 

88 

% Error 
UAI 

Ver. 10.0 

MSC 
Ver. 65C 

1 

3.571x10-5 

(1.406x10-3) 

0.17 

0.13 

0.04 

2 

8.865x10-5 

(3.490x10-3) 

0.14 

0.10 

0.03 

5 

11.34x10-5 

(4.465x10-3) 

0.11 

0.11 

0.07 

10 

11.38x10-5 

(4.482x10-3) 

0.08 

0.08 

0.05 
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TABLE 4: TRIA3 Error in Displacement at Center of Plate 

Aspect Ratio Study (12 x 12 Mesh) 

Stiff Core, Simply-Supported, Honeycomb Plate 
Under Uniform Pressure Load with Transverse Shear Flexibility 


AR 

theoretical w, 
m (in.) 

Cosmic 

88 

% Error 
UAI 

Ver. 10.0 

MSC 
Ver. 65C 

1 

2.422x10-3 

(9.535x10-1) 

0.37 

0.33 

0.24 

2 

5.974x10-3 

(2.352x10-1) 

0.28 

0.24 

0.17 

5 

7.631x10-3 

(3.004x10-1) 

0.21 

0.21 

0.17 

10 

7.660x10-3 

(3.016x10-1) 

0.21 

0.21 

0.17 
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TABLE 5: TRIA3 Error in Displacement at Center of Plate 

Aspect Ratio Study (12 x 12 Mesh) 

Flexible Core, Simply-Supported, Honeycomb Plate 
Under Uniform Pressure Load with Transverse Shear Flexibility 


AR 

theoretical w, 
m (in.) 

Cosmic 

90 

% Error 
UAI 

Ver. 11.1 A 

MSC 
Ver. 66A 

1 

3.102x10-3 

(1.221x10-1) 

0.30 

0.28 

0.20 

2 

7.026x10-3 

(2.766x10-1) 

-0.67 

0.24 

0.18 

5 

8.785x10-3 

(3.459x10-1) 

0.22 

0.22 

0.17 

10 

8.815x10-3 

(3.470x10-1) 

0.17 

0.17 

0.13 
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TABLE 6 


SUMMARY OF TEST RESULTS FOR TRIA3 SHELL ELEMENTS 


Test 

| Elem. Loading 

Element 

Shape 

COSMIC 

90 

UAI 

11.1A 

MSC 

66A 

In 

Plane 

Out of 
Plane 

1. Patch Test 

X 


Irregular 

A 

A 

A 

2. Patch Test 


X 

Irregular 

A 

A 

A 

3. Straight Beam, Extension 

X 


AH 

A 

A 

A 

4. Straight Beam, Bending 

X 


Regular 

F 

F 

F 

5. Straight Beam, Bending 

X 


Irregular 

F 

F 

F 

6. Straight Beam, Bending 


X 

Regular 

B 

B 

B 

7.. Straight Beam, Bending 


X 

Irregular 

B 

B 

B 

8. Straight Beam, Twist 



All 

F 

F 

F 

9. Curved Beam 

X 


Regular 

F 

F 

F 

10. Curved Beam 


X 

Regular 

F 

F 

F 

11. Twisted Beam 

X 

X 

Regular 

C 

C 

D 

12. Rectangular Plate (N=4) 


X 

Regular 

B 

B 

B 

13. Scordelis-Lo Roof (N=4) 

X 

X 

Regular 

D 

D 

D 

14. Spherical Shell (N=8) 

X 

X 

Regular 

A 

A 

A 

15. Thick-Walled Cylinder 

X 


Regular 

A 

A 

A 

(nu=.4999) 







Number of Failed Tests (D's and F's) 


6 

6 ; 

7 


Grading for Shell Element Test Results 


Grade 


Requirement 


A 

B 

C 

D 

F 


2% > Error 
10% > Error > 2% 
20% > Error > 10% 
50% > Error > 20% 
Error > 50% 
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Fig. 1 

Test Problem 



Plate Size: a=1.016 m (40. in.)* b=1.016m (40.in.) 
Boundary Conditions: simply supported on all edges 
Loading: pressure load, p=6895. N/m 2 (1.0 psi) +Z direction 
Thickness: t=0.0508 m (2.0 in.) 

*: Variable in aspect ratio studies 


39 


Fig. 2 

Mesh Geometry 



AR C = a/b e = element aspect ratio 

N x = a/2a e = number of elements in X direction in 1/4 of plate 
N y = b/2b c = number of elements in Y direction in 1/4 of plate 
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Fig. 3a 


Theoretical Solution - Central Displacement 


Central Displacement 


w(x= y,y=0) = J [ 1 + C 5 cosh( M- y) + H yC 6 sinh( V y) 


W 


( j 2+vX-, sin px 

\C S ' C n /J p5 


where. 


Cc = 


1 •- /I 1+ v \ 1 ,, 

[l + p2D(^-- — ) + T « m tanh( ttm 

<Xm S n 


cosh 


c 6 = 


1 


2 cosh am 


mjt b t mt 

a m = ~2~~a > M- = a 
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Fig. 3b 


Theoretical Solution - Homogeneous Plate Parameters 


Homogeneous Plate 


D = 


Et^ 

12(1- v 2 ) 


C 


n 


5 Et 

6 v 



G= 


E 

2(1+ v) 


E = 6.89 x 10 10 N/m 2 (10.0 x 10 6 lb/in 2 ) 
v = 0.33 

t = .0508 m (2.0 in.) 
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Fig. 3c 


Theoretical Solution - Honeycomb Plate Parameters 


Honeycomb Plate 

Eftf (t c +t ) 2 

4(1- V 2 ) 

C n = oo 

C s =t c^c 

E f = 6.89 x 10 10 N/m 2 
(10 x 10 6 lb/in 2 ) 

v = 0.33 

G c = 1.379 x 10 7 N/m 2 (2000. lb/in 2 ) Flexible Honeycomb Plate 
or 

1.517 x 10 8 N/m 2 (22000. lb/in 2 ) Stiff Honeycomb Plate 

t c = .0508 m (2.0 in) 
t f = .254 mm (.01 in.) 

1 = tc + tf 
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Fig. 4 

TRIA3 Error in Displacement at Center of Plate 

Mesh Size Study (Element Aspect Ratio 1 .0) 
Simply-Supported, Homogeneous Plate 
Under Uniform Pressure Load with Transverse Shear Flexibility 



% JO J31U93 JB }U3UI03B|dSja u l JOJJfl 
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Mesh Size, Nx and Ny 



TRIA3 Error in Displacement at Center of Plate 

Mesh Size Study (Element Aspect Ratio 1 .0) 
Simply-Supported, Homogeneous Plate 
Under Uniform Pressure Load without Transverse Shear Flexibility 
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Mesh Size, Nx and Ny 





Fig. 6 

TRIA3 Error in Displacement at Center of Plate 

Mesh Size Study (Element Aspect Ratio 1.0) 

Stiff Core, Simply-Supported, Honeycomb Plate 
Under Uniform Pressure Load with Transverse Shear Flexibility 
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Mesh Size, Nx and Ny 



Fig. 7 

TRIA3 Error in Displacement at Center of Plate 

Mesh Size Study (Element Aspect Ratio 1 .0) 

Flexible Core, Simply-Supported, Honeycomb Plate 
Under Uniform Pressure Load with Transverse Shear Flexibility 
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Mesh Size, Nx and Ny 



Fig. 8 

TRIA3 Error in Displacement at Center of Plate 

Aspect Ratio Study (12 x 12 Mesh) 

Homogeneous, Simply-Supported Plate 
Under Uniform Pressure Load with Transverse Shear Flexibility 
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Element Aspect Ratio 




TRIA3 Error in Displacement at Center of Plate 

Aspect Ratio Study (12 x 12 Mesh) 

Stiff Core, Simply-Supported, Honeycomb Plate 
Under Uniform Pressure Load with Transverse Shear Flexibility 
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Element Aspect Ratio 





Fig. 10 

TRIA3 Error in Displacement at Center of Plate 

Aspect Ratio Study (12 x 12 Mesh) 

Flexible Core, Simply-Supported, Honeycomb Plate 
Under Uniform Pressure Load with Transverse Shear Flexibility 
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Element Aspect Ratio 
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VALIDATION OF THE CQUAD4 ELEMENT FOR VIBRATION AND SHOCK 
ANALYSIS OF THIN LAMINATED COMPOSITE PLATE STRUCTURE 

by 

Douglas E. Lesar 

Ship Structures and Protection Department 
Carderock Division, Naval Surface Warfare Center 
Bethesda, Maryland 20084-5000 U.S.A. 

ABSTRACT 

The CQUAD4 thin plate element implemented in COSMIC NASTRAN is 
capable of modeling thin layered plate and shell structures composed 
of orthotropic lamina. Fiber-reinforced composites are among the 
classes of inhomogeneous and non— isotropic materials which can be 
treated. Although the CQUAD4 has been extensively checked in static 
cases, little validation has been carried out for vibration response 
modeling. This paper documents validation of the CQUAD4 element's 
accuracy for vibration response analysis of thin laminated composite 
plates. 

The lower-order natural frequencies and mode shapes of ten 
glass fiber-reinforced plastic (GFRP) and carbon fiber- reinforced 
plastic (CFRP) plates are computed and compared to published experi- 
mental and numerically-computed data. A range of ply geometries 
including unidirectional, cross-ply, and angle-ply are considered. 
The plates' length-to-thickness ratios all lie in the vicinity of 
100 to 150. The CQUAD4 plate idealizations provide natural frequen- 
cy predictions within ten percent of measured data for all six 
lowest modes of seven of ten plates. For two of the remaining three 
plates, only the fundamental frequency is predicted with an error 
greater than ten percent. Results for the one remaining plate do 
not correlate with published data, possibly because of erroneous 
reporting of its geometry or material properties in the literature. 
To obtain accurate frequency predictions, lamina in-plane elastic 
moduli had to be tuned to reflect each plate's fiber volume 
fraction. 

These results show that the NASTRAN CQUAD4 plate element is 
useful and reasonably accurate for vibration and shock analysis of 
structures composed of thin fiber-reinforced plastic plates. 

INTRODUCTION 

There is strong current Navy interest in exploitation of fiber- 
reinforced plastics as lightweight materials for a wide variety of 
ship structures. These structures must be designed to withstand in- 
service loads of quasi-static, transient dynamic, and steady-state 
dynamic nature. For many ship structures, transient shock is a pri- 
mary load. For these and others, steady-state vibration response 


51 


impacts the ship's acoustic signature. In both cases, the struc— 
ture's modal properties (natural frequencies, mode shapes, and modal 
loss factors) are key parameters governing its transient and steady- 
state response. In most cases, the ship structures being designed 
are complex enough so that numerical (finite element) methods must 
be employed to obtain realistic modal property estimations. 

The NASTRAN finite element code is one of the Navy's premier 
tools for steady-state vibration response analysis of ship and 
submarine structures. Undamped natural mode analysis and forced 
vibration response analysis with hysteretic damping can be performed 
by NASTRAN, as well as modal frequency and loss factor analysis for 
structures with viscoelastic damping materials (Ref. 1) . NASTRAN is 
also a key component of the NASHUA suite of codes for performing 
radiated noise and acoustic scattering analysis of vibrating submer- 
ged structures (Ref. 2,3). 

The COSMIC NASTRAN CQUAD4 element is designed to model aniso- 
tropic layered plates as well as homogeneous and isotropic plates. 
The theory and assumptions behind the CQUAD4 have been informally 
documented (Ref. 4). The CQUAD4 has been found to be more accurate 
for a given finite element grid than its predecessor, the CQUAD2 
element, for prediction of low frequency eigenmodes of thin-walled 
cylindrical shells composed of isotropic materials (Ref. 5) . To the 
author's knowledge, no comparable study of the accuracy of the 
CQUAD4 formulation for anisotropic plate or shell vibration yet 
exists, particularly for plates and shells composed of layers of 
fiber-reinforced plastic lamina. 

This paper summarizes an investigation of the accuracy of the 
NASTRAN CQUAD4 membrane and plate bending element for vibration 
analysis of structures composed of thin fiber-reinforced composite 
plates. This is accomplished by comparisons of NASTRAN-computed 
undamped natural frequencies of ten GFRP and CFRP plates with 
published experimental data and other numerical predictions. 

NOTATION 

E X1 Lamina extensional modulus in 

direction parallel with fibers 

E 22 Lamina extensional modulus in 

direction transverse to fibers 

Through-thickness extensional 
modulus 

G 12 In-plane lamina shear modulus 

G 13 Transverse lamina shear modulus for 

out-of-plane shearing of a fiber 
cross section 
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G 2 3 

Transverse lamina shear modulus for 
in-plane shearing of a fiber cross 
section 

L 

Plate side length 

t 

Plate thickness 

V f 

Fiber volume fraction 

V 12 

In-plane shear strain 

v V 

r 13' '23 

Out-of-plane shear strains 

e x , € 2 

In-plane extensional strains 

e 

Fiber orientation angle with 
respect to one plate side 

^12' ^21 

In-plane lamina Poisson's ratios 

^ 23 

Transverse lamina Poisson's ratio 

c 

Average mass density 

*~1' **2 

In-plane extensional stresses 

/yo 
' 12 

In-plane shear stress 

13' 23 

Out-of-plane shear stresses 


ABBREVIATIONS 

avg. 

Average 

CFRP 

Carbon fiber-reinforced plastic 

CLT 

Classical lamination theory 

consis. 

Consistent 

CPT 

Classical plate theory 

DOF 

Degree (s) -of-freedom 

FEA 

Finite Element Analysis 

GPa 

Gigapascals 

GFRP 

Glass fiber-reinforced plastic 

Hz 

Hertz 

kg 

Kilograms 
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m 

Meters 

mm 

Millimeters 

params • 

Parameters 

psi 

Pounds (force) per square inch 

ref • 

Reference (s) 

RMS 

Root-mean-square 

sym 

Midplane-symmetric 


BACKGROUND 

Anisotropic plate analysis is more difficult and involves more 
variables and parameters than isotropic plate analysis. The impor- 
tance of accounting for transverse shear flexibility in relatively 
thin composite laminates is discussed in the next section. The ways 
by which the CQUAD4 addresses these difficulties are also briefly 
described. 

TRANSVERSE SHEAR DEFORMATION ISSUES 

In composite plate mechanics, the counterpart of the well-known 
Kirchhoff (Classical plate) theory (CPT) for isotropic plates is the 
so-called "Classical Lamination Theory" (CLT) . The two theories 
invoke the same kinematic assumptions regarding the deformation of 
the plate with respect to its middle surface; that is, sections 
originally planar and perpendicular to the middle surface remain 
planar and perpendicular in the deformed state. The mathematical 
development of CLT occupies much of the text by Jones (Ref. 6) . The 
reader should consult this (or some other) text for detailed 
exposition of CLT assumptions and derivation of the CLT equations. 

The accuracy of the Kirchhoff kinematic assumption in isotropic 
piste mechanics degrades when the plate thickness becomes signifi- 
cant compared to its span length. For isotropic metallic plates, 
transverse shear-stiffness-governed transverse deflection becomes 
significant relative to flexural deflection when the span length-to- 
i-h.ic-kness ratio (L/t) is sufficiently small. An idea of required 
smallness can be gained from the discussion of "corrected" plate 
flexural waves in Chapter II, Section 3b of Cremer, Heckl, and 
Ungar (Ref. 7) . They show, for a uniformly thick plate composed of 
isotropic material, that transverse shear effects decrease flex- 
ural wavespeed by ten percent when the flexural wavelength is about 
six times the plate thickness. For one-half wavelength over a plate 
span length, this limitation translates to an L/t ratio of 3. 

In the case of high modulus composite plates, an analogous 
limitation of the adequacy of CLT is encountered at larger L/t. 

This more stringent limit arises because the ratio of effective 
laminate extensional elastic modulus to shear modulus is a key fac- 
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tor (other than section geometry) governing the magnitude of shear 
deformation relative to flexure. A rough idea of limiting L/t 
ratios for the adequacy of CLT for composite plates follows from 
insertion of some "ballpark" ratios of effective laminate extension- 
al and transverse shear moduli for GFRP and CFRP into the approxi- 
mate expression for shear-corrected plate flexural wavespeed found 
on page 115 of Reference 7. Using E.. / G 2 _ = 13 for GFRP and 125 
for CFRP, ten percent differences in wavespeed arise for L/t of 6.75 
and 20.9, respectively. (A comparable modulus ratio . for isotropic 
materials is 2.6). These ad hoc assessments are qualitatively corro- 
borated by static examples found in section 6.5 of Reference 6. In_ 
the problem of cylindrical bending of a CFRP strip with E 1 - / G - 
125, maximum static deflection predicted by. CLT is twenty percent 
smaller than the true shear-corrected solution at L/t = 20. 

The low frequency composite plate vibration literature is domi- 
nated by evaluation of methods for account of transverse (interlami- 
nar) shear in prediction of natural vibration frequencies and modal 
deflections. The inadequacies of CLT even for fundamental plate 
frequencies are repeatedly demonstrated in the literature for. L/t s 
of 5 or 10. Many approaches have been developed to provide finite 
element-based plate formulations to handle through-thickness stress 
fields for arbitrary L/t. In these complex formulations, transverse 
normal stresses are no longer assumed to be zero, transverse shear 
stresses are constrained to be zero on upper and lower laminate 
surfaces, and shear stress continuity between laminae is maintained. 
Ag a result, sections perpendicular to the plate midplane rotate, 
with respect to the midplane and warp out of a planar configuration. 

The plate finite elements simulating through-thickness stress 
fields in laminates are, obviously, quite sophisticated. Tables III 
and VI of Mallikarjuna and Kant (Ref. 8) provide a good flavor for 
the performance of some higher-order approaches for reckoning with 
the effect of transverse shear and normal strains in plate vibration 
frequency prediction. Table III shows that CLT provides reasonable 
fundamental frequency predictions for a simply- supported CFRP angle- 
ply plate for L/t at and above 20. Even for such highly anisotropic 
plates, the complications of very complex high-order transverse 
shear theories, necessary for laminate strength and structural 
integrity problems, are not justified in vibration problems unless 
L/t is less than about 20, or if short-wavelength (high frequency) 
vibration modes are of interest. 

The NASTRAN CQUAD4 membrane and bending element embodies the 
assumptions of CLT, but contains first-order corrections for trans- 
verse shear flexibility. When the CQUAD4 is used to model homogen- 
eous plates, transverse shear strains vary linearly with the thick- 
ness coordinate, and zero shear stress and strain boundary condi- 
tions on upper and lower plate surfaces are not satisfied. The 
shear energy implied by the linear distribution is corrected by a 
multiplicative constant to produce the energy implied by the true 
quadratic distribution. As a result of these assumptions, sections 
perpendicular to the plate midplane are allowed to rotate with 
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respect to the midplane, but remain straight when the plate is in 
the deflected state. This ad hoc transverse shear correction ex- 
tends the element's range of validity to thicker homogeneous plates, 
but is not acceptable for inhomogeneous layered plates. 

When modeling inhomogeneous plates composed of orthotropic 
layers, a quadratic transverse shear strain distribution is assumed 
in each layer of the CQUAD4 element. Interlaminar shear strains are 
matched at lamina interfaces, and zero shear stress boundary condi- 
tions are enforced on the upper- and lower-most lamina surfaces. 
However, the through-thickness normal stress is assumed to be zero, 
and complete consistency between the strain-displacement equations 
for in-plane direct strain and transverse shear strain is not main- 
tained. The CQUAD4 is thus seen to overcome the limitations of CLT 
for laminated plates, but does not represent all aspects of the 
kinematics of three-dimensional continua taking place in relatively 
thick laminates. 

The CQUAD4 element is discussed in more depth in the following 
section. 

THE CQUAD4 ELEMENT 

The CQUAD4 is a four-noded planar element possessing membrane, 
flexural, and transverse shear stiffness. In the case of layered 
plates, individual laminae are not modeled explicitly; rather, 
equivalent stiffness matrices for the plate as a whole are defined. 
Each lamina is assumed to be in a state of plane stress, and the 
laminae are presumed to be perfectly bonded by infinitesimally thin 
non-shear-deformable layers. Each lamina is assumed to be specially 
orthotropic, with six independent elastic moduli when through- 
thickness direct stresses and strains are ignored. Any alignment of 
lamina fiber axis with respect to the local element coordinate 
system can be accomodated. Hence, any layup or stacking sequence 
can be handled (unidirectional, cross-ply, regular or irregular 
angle-ply) . Layups unsymmetric with respect to plate midplane can 
also be modeled, as membrane-bending coupling is accounted for when 
it occurs. Each lamina may also be composed of a different ortho- 
tropic material, if desired. 

The element's stiffness matrix terms for determining in-plane 
displacements and flexural rotations as a function of imposed forces 
and moments arise from CLT assumptions. The CQUAD4's force-versus 
strain equations for membrane, flexural, and membrane-flexure coup- 
ling (Ref. 4) are identical to those developed in sections 2.1 
through 2.6 of Reference 6. The kinematic assumptions regarding in- 
plane and flexural strains and displacements follow classical 
assumptions and require no explanatory remarks here. Reference 4 
documents the force-versus strain matrix terms for transverse shear 
strains. These are based on overall element equilibrium, continuity 
of transverse shear between adjacent laminae, and satisfaction of 
zero shear strain and stress boundary conditions on the upper- and 
lower-most laminate facings. As mentioned earlier, the strain- 
versus displacement matrix is based on the assumption that through- 
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thickness planar sections rotate with respect to the plate midsur- 
face and also distort out of planes originally perpendicular to the 
midsurface. 

The CQUAD4 is an isoparametric element, whose displacement 
fields are interpolated through space by linear variations of in- 
plane and transverse displacements (and rotations about midsurface) 
between grid points. The associated in-plane strains are constant 
between grid points but vary linearly with the thickness coordinate. 
However, transverse shear strain varies quadratically through each 
lamina thickness. 

Those who use the CQUAD4 to model fiber-reinforced composite 
plates must realize that the input elastic moduli are effective 
moduli for a particular fiber-matrix combination . There are many 
different fibers in use (glass, carbon, kevlar, and boron are exam- 
ples) and many resins or matrix materials, each of which has their 
own unique elastic moduli. Although matrix resins are usually con- 
sidered to be isotropic, fibers have distinct extensional, trans 
verse, and shear moduli. The effective in-plane moduli of an ortho- 
tropic continuum defined to be equivalent to the actual inhomogen 
eous fiber and resin system are (sometimes nonlinear) functions of 
the extensional, transverse, and shear moduli of the fiber, the 
extensional and shear moduli of the resin, and the Poisson s ratios 
of the fiber and resin. The fraction of the lamina volume occupied 
by fiber material is an important variable defining the magnitude of 
effective in— plane lamina moduli. Some strength— of— materials rela- 
tionships defining effective lamina in-plane moduli as a function of 
constituent matrix and fiber moduli and fiber volume fraction are 
developed in sections 3.1 and 3.2 of reference 6. Reference 9 pro- 
vides a handy tabulation of effective lamina elastic moduli under 
the assumption of a transversely isotropic lamina. Chamis provides 
formulas for effective out-of-plane shear moduli as well as the more 
commonly reported in-plane extension and shear moduli. 

APPROACH 

Specifics of the present vibration modeling study are now 
described. 

PLATE VIBRATION SPECIMENS 

It was desired that the validity of the CQUAD4 be proven by 
comparison to numerical results obtained independently by other 
researchers, and to experimental data, if possible. Lin, Ni, and 
Adams (Ref. 10,11) and Xiao, Lin, and Ju (Ref. 12) have published 
experimental data and finite element computations for the lowest six 
vibration modes of square CFRP and GFRP plates for free, uncon- 
strained boundary conditions. They measured and computed both natu- 
ral frequencies and damping loss factors, which makes their work 
almost uniquely complete and thorough. Reference 11 contains data 
for nine plates repeating that for four plates treated in reference 
10. Reference 12 revisits three previously examined plates with a 
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more refined finite element formulation allowing plane sections to 
rotate and warp, with consistent correction of all strain-displace- 
ment relations for these kinematic conditions. This paper contains 
data for one additional new plate, making a total of ten unique 
plates in all (four CFRP and six GFRP) . Refinement of their finite 
element method (Ref. 12) improved correlation of their predictions 
with measurements. 

The geometric and material parameters of the ten plates studied 
by Lin et. al. and Xiao et. al. are listed (in SI units) in Table 1. 
The parameters given in reference 12 for plates 770 and 772 are 
inconsistent with those reported in reference 11. The NASTRAN study 
confirmed that the side lengths, thicknesses, average mass 
densities, and fiber volume fractions given for plates 770 and 772 
in reference 11 are the correct values. Further correlation of 
plate parameters, frequencies, and mode shapes between references 
10-12 revealed accidental reversal of mode shape plot labels in ref- 
erence 11. In addition, two unidirectional GFRP plates of different 
size are reported in references 11 and 12 with the identification 
number 761. (They are herein distinguished from the other as 761L 
and 761X) . These discrepancies initially caused much confusion, but 
the author is confident that the data in Table 1 is correct. 

Some features of these plate specimens are notable. All of 
them are square, and have L/t ratios in the vicinity of 100 to 150. 
Even though such L/t would seem to be in the range of applicability 
of CLT, Xiao et. al. show that CLT, which totally ignores transverse 
shear deformations, overestimates natural frequencies by factors as 
high as sixteen percent over a theory with first-order shear correc- 
tion. All of the laminates listed in Table 1 are symmetric about 
the plate midplane, eliminating flexure/extension coupling effects. 
Five plates (GFRP specimens 734, 761(L), 761(X) and CFRP specimens 
762 and 764) have "specially orthotropic" lamina (all fibers aligned 
with the plate sides). For these plates, there is no coupling 
between in-plane extension and in-plane shearing, so their vibration 
mode shapes have nodal lines more or less parallel with the plate 
edges. The other five plates (765, 769, 771, 770, 772) have at 
least some plies with fibers angled relative to plate edges, pro- 
viding more complex mode shapes. 

All plates were tested with "free" edges (supported by soft 
foam rubber strips) and numerically analyzed with zero-constraint 
boundary conditions. (Although many practical design problems of 
interest would involve plate structures with boundary constraints, 
this feature eliminates uncertainties about which boundary degrees- 
of-freedom (DOF) to constrain to obtain "simply supported" edges) . 

An iterative technique was used to obtain natural frequencies and 
modal damping loss factors. First, the specimens were excited into 
steady-state vibration by an electrodynamic shaker, and approximate 
resonance frequencies were determined. Nodal lines for each excited 
natural mode were located by the classical Chladni sand pattern 
technique. Then, for each mode of interest, locations of the 
support strips were then adjusted to align with nodal lines. A 
transient excitation technique was then used to obtain more precise 
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estimates of resonance frequencies and modal damping loss factors. 
LAMINA MATERIAL PROPERTIES 


Effective in-plane lamina elastic moduli are reported in Lin 
et. al. and Xiao et. al. These are listed in Table 2, in both SI 
and English units. These moduli pertain to a lamina microstructure 
where half of the total ply volume is occupied by fibers and half by 
the matrix. For each plate analyzed here, the moduli input to 
NASTRAN must be adjusted up or down from these nominal values, 
according to the measured fiber volume fraction (V f ) for each plate. 
In general, the moduli vary nonlinearly with V_, but experimentally- 
verified semi-empirical equations are available for determining V~- 
adjusted moduli. References 10-12 reported only the nominal moduli 
in Table 2, but provided a literature source (Ref. 13) for accom- 
plishing the adjustments. 

For orthotropic laminae in a state of plane stress (in the 1-2 
plane) , the only independent engineering moduli are E ll' - 22 ' ? 12 ' 
and V 12 - The additional Poisson's ratio V 21 must satisfy the 
relationship: 

V 12 / E 11 = ^21 / E 22 

The NASTRAN CQUAD4 element enforces this constraint on V (Ref. 

4) . With these moduli, all in-plane strain and stress components 
are defined by the CLT. However, the NASTRAN CQUAD4 element also 
requires, as input, nonzero transverse shear moduli G-_ and G_ 3 , 
associated with out-of-plane shearing, which are not Specified or 
required in CLT. Fortunately, fiber-reinforced lamina may often be 
assumed to be "transversely isotropic" in analysis of composite 
structures, as implied by reference 9. This means that if the 
lamina lies in the 1-2 plane, and the fibers are aligned in the 1 - 
direction, then the transverse shear modulus G__ is equal to the 
in-plane shear modulus G. _ . The transverse shear modulus associated 
with shearing of the matrix "around" the fibers and distortion of 
the fiber in a plane perpendicular to its axis (<* 23 ) remains to be 
determined. To clarify these matters, the stress and strain compo- 
nents for a transversely isotropic lamina with zero direct stress 
normal to the lamina are illustrated in Figure 1. 


Lin et. al. and Xiao et. al. did not report either of the 
transverse shear moduli used in their analyses. Educated guesses 
had to be made for both the CFRP and GFRP G . Reference 8 provides 
nondimens ional elastic moduli for a CFRP— liRe material which are 
closely satisfied by the in-plane parameters for the CFRP in Table 
2. According to reference 8 , G _ 3 should equal 0.2E__. This implied 
value of G-- was used in the present CFRP plate analyses, based on 
fiber volume fraction-adjusted E_ 2 . For the glass fiber lamina, G 23 
= 0.6G was assumed, based on in— house experience with such mater- 
ials. ^?he assumed out— of— plane moduli are indicated in the "notes" 
section of Table 2. G 23 could also be estimated by methods reported 

in reference 9 . 
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Table 3 contains lamina elastic moduli for each of the ten 
plates, after adjustment of the nominal moduli in Table 2 for fiber 
volume fraction. Lin et. al. did not correct the in-plane Poisson's 
ratio for fiber volume fraction, even though Ni and Adams show that 
it decreases appreciably for increasing V f . Lin et. al. claim that 
accounting for this decrease had no significant effect on their 
calculated natural vibration frequencies. We match their assumption 
by keeping constant at 0.3 for all ten plates. In general, 

V 12 should be adjusted for V f . 

ASSUMPTIONS REGARDING LAMINATE BEHAVIOR 

The kinematic assumptions of the CQUAD4 element have been dis- 
cussed previously; namely, that planar sections perpendicular to the 
plate midplane can rotate and warp relative to the midplane. NAS- 
TRAN-computed natural frequencies for CQUAD4 plate idealizations 
will, subsequently, be compared to some CLT predictions (Ref. 12) 
and to the FEA predictions of Lin et. al. (Ref. 11) and Xiao et. al. 
(Ref. 12) These three approaches differ in accuracy, and it is 
important to understand how the NASTRAN predictions should compare 
with them. Kinematic assumptions are discussed first. 

As discussed previously, CLT totally ignores any transverse 
shear deformation effect, and will always predict lower-order 
vibration frequencies which are too high for laminates with low L/t 
and high extensional-to-shear modulus ratios. In terms of through- 
thickness lamina kinematics, CLT prescribes linear variations of 
direct and in-plane shear strains, and zero transverse shear 
strains. The kinematic assumptions of Lin et. al. (Ref. 11) and 
within the CQUAD4 (Ref. 4) are identical, and imply linear through- 
thickness variations of direct and in-plane shear strains and 
quadratic through— thickness variations of transverse shear strains, 
with zero shear stress boundary conditions on the upper- and lower- 
most lamina facings satisfied. However, this thick plate-type 
theory is only an approximation; consistent correction of all 
strain-displacement relationships when cross-section warping is 
allowed requires cubic through-thickness variations of direct and 
in-plane shear strains. Xiao et. al. (Ref. 12) include this consi- 
derable complication in their plate-type FEA formulation, which can 
be understood as a special case of three-dimensional elasticity. 

From an understanding of kinematic assumptions alone, CQUAD4 
composite plate idealizations should provide natural frequency 
estimates that are (1) , more accurate than CLT, (2) as accurate as 
the Lin et. al. predictions, and (3) less accurate than the Xiao et. 
al. predictions. However, other factors will influence the com- 
parisons between FEA-predicted frequencies; particularly, the 
polynomial form of the interpolation functions expressing element 
displacement fields in terms of grid point displacements, and the 
mass matrix formulation (consistent or lumped) . The CQUAD4 utilizes 
linear interpolation functions with respect to the element's four 
grid points. Lin et. al. do not specify the interpolation func- 
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tions used in their M 8-node 40 -degree-of-freedom" elements, although 
they can be guessed to be of at least guadratic order. 

NASTRAN ANALYSES 

The results of the NASTRAN computations are now compared to 
the calculations and experimental data in references 10-12. 

FEATURES OF THE NASTRAN PLATE IDEALIZATIONS 

Lin et. al. and Xiao et. al. employed a 6-by-6 mesh of 8-noded 
isoparametric rectangular elements in all of their plate FEA ideali- 
zations. Each element possessed 40 DOF, (5 per grid point), imply- 
ing that all rotations about axes normal to the plate surface were 
constrained. As in the CQUAD4 , they applied a numerical condition- 
ing factor to transverse shear stiffness terms to eliminate excess- 
ive shear stiffness ("shear locking") ; a consequence of numerical 
integration of element stiffness. The NASTRAN CQUAD4 elements also 
have five DOF per node, with normal axis rotation constraints app- 
lied, and are also conditioned to avoid shear locking. 

The CQUAD4 differs from the Lin and Xiao et. al. elements in 
one important way; they are 4-noded isoparametric quadrilaterals and 
thus have linear interpolation of in-plane displacements between 
grid points instead of quadratic interpolation. Thus, although the 
Lin et. al. element is kinematically similar to the CQUAD4 as far as 
through-thickness shear effects are concerned, the elements' assumed 
in— plane displacement fields differ as a function of the plate's in- 
plane dimensions. The Xiao et. al. element provides for a more com- 
plex through— thickness displacement and strain distribution than the 
Lin et . al. element and the CQUAD4, as discussed earlier. 

The NASTRAN CQUAD4 meshes used here consist of a 12-by-12 grid 
of elements, with 169 grid points and 845 unconstrained DOF. These 
idealizations are roughly comparable to the Lin and Xiao models in 
modal displacement field interpolation quality. Three of the plates 
were initially modeled with 6— by— 6 meshes with 49 grid points and 
245 DOF, but these idealizations did not provide sufficient mode 
shape resolution to be acceptable. The Lin and Xiao et. al. mesh is 
compared to the NASTRAN mesh in Figure 2. 

The edges of the NASTRAN plate idealizations were unconstrain- 
ed, and all rotations about axes normal to the plate surface were 
suppressed. A SUPORT input record imposing fictitious constraints 
on all five remaining DOF on one grid point was used to provide 
zero— frequency rigid body modes. The FEER eigenmode extraction^ 
method was used in NASTRAN Rigid Format 3 (normal modes analysis) , 
with the lowest twenty modes requested. The Inverse Power method 
was employed in some trial runs; it provided about half as many 
frequencies as FEER but at more than two times greater run time and 
cost. All computations were performed by RPK COSMIC NASTRAN, 1990 
release, installed on the DTRC CRAY X-MP supercomputer in a COS 
operating system environment. Typically, 22 to 24 modes were ex- 
tracted by FEER in 23 to 24 CP seconds. 
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Lin and Xiao do not mention whether they used lumped or con- 
sistent mass matrix formulations. Both options were considered in 
the NASTRAN study. 

NASTRAN MODAL ANALYSIS RESULTS 

Initial NASTRAN analyses utilized nominal elastic moduli based 
on fifty percent fiber volume fraction. Although predicted frequen- 
cies were in the general vicinity of measured values, they did not 
compare with CLT results, the Lin et. al. results, or the Xiao et. 
al. results in the expected way. That is, the CQUAD4 did not neces- 
sarily appear more accurate than CLT and less accurate than the Xiao 
et. al. approach. It suffices to say that fiber volume fraction 
strongly influences effective moduli and must be accounted for to 
obtain credible composite plate modal predictions. 

Relative NASTRAN frequency predictions are compared with CLT, 
Lin's, and Xiao's results in Table 4. There is, in general, good 
correlation with measurements, and NASTRAN 's frequency error trends 
are comparable to those in the Lin et. al. analyses (especially 
those for a consistent mass matrix). The fact that effective lamina 
elastic moduli must be corrected for fiber volume fraction to obtain 
credible natural frequency predictions for laminated composite 
plates is emphasized in Table 5, where root-mean-square (RMS) values 
of frequency prediction percentage error are tabulated for the six 
modes of each plate. RMS error is seen to be significantly reduced 
in most cases when fiber volume fraction-corrected lamina moduli are 
employed. The major exception is plate 765, which still suffers 
from some large unknown systematic error making predicted frequen- 
cies far too low. Results for plate 770 were not much changed since 
actual V f was already close to one-half. RMS errors for plates 764 
3i"id 771 remain at ten percent and above. However, the major part of 
these high errors involves their fundamental modes, which can be 
easily impacted by test article boundary constraint. (The Lin et. 
al. analysis also predicted overly high frequencies for these 
modes) . 

Absolute measured and predicted natural frequencies are summar- 
ized in Table 6, which lists Lin et. al. measurements and computa- 
tions and the NASTRAN CQUAD4 idealization results for lumped and 
consistent mass. Mode— by— mode and average percentage frequency 
prediction errors are listed in Table 7 along with mode shape de- 
scriptions . The NASTRAN lumped mass model gives the lowest average 
error in seven of ten cases, and the NASTRAN consistent mass model 
is best in two of the three remaining cases. (Interestingly, the 
consistent mass model errors roughly parallel those of the Lin 
model). Lin et. al. calculations are "best" only for plate 765, for 
which the author believes a parameter was misdocumented in reference 
11. In six of the nine plates other than 765, NASTRAN predicts 
frequencies that are all at or within ten percent of measured. In 
the other three, errors are greater than ten percent only for the 
first and/or second modes, and these roughly parallel the errors in 
the Lin et. al. analyses. 
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Plots of chosen NASTRAN-computed eigenmodes are presented in 
Figures 3-8. Plate 761X does not appear since its mode shapes are 
identical to plate 761L, and plates 764, 770, and 772 are omitted 
since their mode shapes are very similar to those of plates 73 4, 

769 and 771. The nodal patterns can be confirmed as being identi- 
cal' to those measured. (Figures 3-8 should be compared to Tables 8, 

7 10 9, 11, and 3, respectively, of reference 11). Finer details 

of some NASTRAN modes, particularly the veering of nodal lines away 
from each other in plates with angled plies, are more easily seen in 
colored graphics terminal displays. Only the plates with specially 
orthotropic ply layups (734, 761L, 761X, 762, 764) exhibit eigen- 
modes with nodal lines more or less parallel to the plate sides. 

Most of these modes are essentially beam-like flexural modes with or 
without phase changes at one symmetry plane. In contrast, plates 
with angled plies possess a larger number of more complex plate 
flexural modes. 

The natural frequency results obtained via NASTRAN CQUAD4 plate 
element idealizations in these simple composite plate vibration 
problems are judged to be acceptably accurate for engineering pur- 
poses. The element performs as well as alternative formulations of 
similar accuracy (the Lin et. al. element) for both GFRP and highly 
anisotropic CFRP plates for a variety of ply geometries. Potential 
users of the element must be cautioned that this validation effort 
concerned plates with L/t ratios in the vicinity of 100 to 150. For 
such L/t, modeling of laminate transverse shear stiffness helps to 
eke out a few percent in low-order natural frequency accuracy , but 
is not absolutely essential to obtain rough-cut results. A more 
critical test of the CQUAD4 for modeling highly anisotropic lami- 
nates would have to concern plates with L/t lower than, say, about 

50. 


It should be mentioned that Lin et. al. and Xiao et. al. also 
measured and computed the specific damping capacities (2 jt times the 
modal loss factor) of each mode. Their FEA program was capable of 
model inq orthotropic lamina damping properties. No attempt was made 
to predict modal damping factors in this effort, as NASTRAN is curr 
ently restricted to the modeling of isotropic material damping. 

SUMMARY 


The performance of the NASTRAN CQUAD4 membrane and plate 
element in analysis of undamped natural vibration modes of thin 
fiber-reinforced composite plates has been evaluated. The element 
provides natural frequency estimates that are comparable in accura- 
cy to alternative formulations, and, in most cases, deviate by less 
than ten percent from experimentally measured frequencies. The 
predictions lie within roughly equal accuracy bounds for the two 
material types treated (GFRP and CFRP) , and for the ply layups 
considered (unidirectional, cross-ply, angle-ply). Effective 
elastic lamina moduli had to be adjusted for measured fiber volume 
fraction to attain this level of accuracy; nominal moduli at fifty 
percent volume fraction gave significantly inferior frequency 
estimates. The lumped mass option provided more accurate frequen- 
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cies than the consistent mass option. 

This evaluation concerned only plates with L/t ratios on the 
order of 100 to 150. Since the CQUAD4 utilizes first-order correc- 
tions for transverse laminate shear stiffness, the element should 
provide useful frequency estimates for plate-like structures with 
lower L/t. For plates with L/t below 20, consideration should be 
given to idealizing with 3-D solid elements. 

Based on the observation that natural frequencies and mode 
shapes are predicted with acceptable engineering accuracy, it is 
concluded that the CQUAD4 should be a useful and accurate element 
for transient shock and steady-state vibration analysis of Naval 
ship structures. 
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Table 1. Geometric and material parameters 
of ten fiber -reinforced plastic 
plates tested and analyzed by 
Lin et. al. and Xiao et. al. 


1 

Plate | Plate 
Ident. | Material 
Number | 

i 

j Plate j 
| Length | 
j L j 

| (ran) | 

l 

Plate | 
Thickness | 
t | 

(mm) j 

Fiber 

Volume 

Fraction 

v f 

Mass 

Density 

<?3 
(kfl/w ) 

Number 

of 

Plies 

| Ply 

| Layup 

| (degrees; see 
| Figure 2) 

| Source 
1 

| Data 
|(ref. no.) 
1 

734 j 

GFRP 

j 227.0 j 

2.05 j 

0.451 

1813.9 

8 

| 10/90/0/90] 

i 

i 

10,11,12 

761 L | 

GFRP 

| 182.75 j 

1.64 | 

0.568 

1971.0 

8 

1 [0] 

I 

11 

761X j 

GFRP 

j 249.0 j 

2.28 j 

0.530 

1924.7 

8 

j CO] 

I 

12 

765 | 

GFRP 

| 230.5 j 

1.45 j 

0.607 

2023.6 

8 

1 C45/-45/45/-45] 

I 

11 

769 

GFRP 

| 224.2 j 

1.37 j 

0.621 

2041.7 

8 

| [0/90/45/-45] 

l 

11 

771 j 

GFRP 

| 204.6 | 

2.11 | 

0.592 

2003.5 

12 

j [(0/-60/60) ] 

I 

10,11 

762 j 

CFRP 

j 178.0 j 

1.58 | 

0.516 

1566.0 

8 

| CO] 

I 

11 

764 j 

CFRP 

1 234.5 j 

2.12 | 

0.342 

1446.2 

8 

j [0/90/0/90] 

I 

10,11 

770 

CFRP 

j 215.0 j 

1.62 j 

0.494 

1551.4 

8 

j 10/90/45/-45] 

i 

11,12 

772 | 

CFRP 

| 215.6 | 

2.02 | 

0.618 

1636.4 

12 

j C(0/-60/60)^] 

i 

10,11,12 


NOTE: All laminates are symmetric about plate midplane. 


Table 2. Nominal in-plane effective elastic 
moduli of GFRP and CFRP laminae at 
fifty percent fiber volume fraction 


Material 


Fiber | Resin | 
type | type | 


22 


12 




i 

1 

| (Gpa) 

1 

<ps £ 

I 

(GPa) 

I 

<ps £ ) 

l 

(GPa) 

I 

| 



i 

1 

I 

l 

i 

1 

/10 ) 

I 


I 

/10 ) 

l 


i 

/10 ) 1 


GFRP 

| "Glass” 

j DX210 

| 37.78 

1 

5.48 

i 

10.90 

i 

1.58 

1 

I 

4.91 

i 

0.71 | 

0.3 


I 

| epoxy 

l 

1 


i 


i 


I 


i 

l 


CFRP 

| HM-S 

j DX210 

| 172.7 

1 

25.0 

i 

7.20 

I 

1.04 

I 

3.76 

I 

0.55 | 

0.3 


1 

| epoxy 

1 

1 


I 


I 


l 


i 

1 


DX210 

1 

I 

| 3.21 

1 

0.47 

I 

3.21 

I 

0.47 

l 

1.20 

i 

0.17 | 

0.34 

epoxy 

1 

I 

1 

1 


I 


I 


l 


l 

1 



NOTE: 


23 

G 23 


0.6(G ) = 2.94 GPa is assumed for GFRP f 

12 

0.2(E ) = 1.44 GPa is assumed for CFRP. 
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Table 3. In-plane and out-of-plane effective elastic 
moduli for GFRP and CFRP laminae, adjusted 
for fiber volume fraction of each plate 



Plate 

Plate 

Fiber 

E 

E , 

l 

i 


G 


Ident. 

Material 

Volume 

1, 

22' 

E„ 

i 

l 

12' 

G 

23 


Number 


Fraction 

(GPa) 

33 

I 

13 

(GPa) 




v f 


(GPa) 

I 

(GPa) 



734 

GFRP 

0.451 



... 

34.4 

9.7 

l 

4.3 

2.6 


761 L 

GFRP 

0.568 

42.5 

13.0 

l 

5.7 

3.4 


761 X 

GFRP 

0.530 

39.9 

12.0 

l 

5.2 

3.1 


765 

GFRP 

0.607 

45.2 

14.2 

l 

6.5 

3.9 


769 

GFRP 

0.621 

46.2 

14.4 

l 

6.7 

4.0 


771 

GFRP 

0.592 

44.1 

14.0 

l 

6.2 

3.7 


762 

CFRP 

0.516 

178.0 

7.4 

l 

3.9 

1.5 


764 

CFRP 

0.342 

119.0 

5.6 

l 

2.6 

1.1 


770 

CFRP 

0.494 

171.0 

7.1 

l 

3.7 | 

1.4 | 


772 

CFRP 

0.618 

213.0 

8.7 

l 

5.0 | 

1.7 | 

1 


NOTE: V 12 - 0.30 is assumed for all laminates, with no adjustment. 


Table 4. Comparison of natural frequencies from final 

NASTRAN C QUAD 4 idealizations and predictions of 
Lin and Xiao et. al. to measured data 


Natural Frequency Ratios, (computed / measured) 


Plate 

Nunber 

and 

pa rams . 


734 
GFRP 
8 plies 

(0/90/0/ 

90]sym 

L/t = 
110.7 


Mode 

Number 


1 

2 

3 

4 

5 

6 


CLT 

ref. 12 
1.20 
1.09 
1.20 
1.11 
1.12 
1.11 


Lin et. 
al. ref. 
11,12 

1.07 

1.00 

1.03 

1.05 

1.04 

1.06 


Xiao 
et. al. 
ref. 12 

0.98 

1.00 

1.00 

1.03 

1.01 

1.01 


NASTRAN 

12-by-12 

lumped 

1.03 

0.94 

0.98 

0.99 

1.00 

0.98 


NASTRAN 

12-by-12 
cons i s . 

1.05 

0.97 

1.01 

1.03 

1.03 

1.04 
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Table 4. (Continued) 


Plate | Mode | Natural Frequency Ratios, (computed / measured) 


and | 

pa rams . | 


1 

1 

1 

1 

CLT 

ref. 12 

Lin 

et. al. 
ref. 11 

l 

l 

1 

Xiao 
et. al. 
ref. 12 

NASTRAN 

12-by-12 

lumped 

NASTRAN 

12-by-12 

consis. 

761 L | 

1 

1 

— 

1.13 

1 

— 

1.11 

1.13 

GFRP j 

8 plies | 

2 

1 

1 

— 

1.00 

l 

l 

| 

— 

0.97 

1.00 

[0] sym | 

3 

1 

1 

— 

1.05 

1 

1 

1 

— 

1.02 

1.05 

L/t = | 

4 

i 

l 

— 

1.00 

1 

l 

— 

0.94 

0.97 

111.4 j 

5 

l 

1 

— 

1.04 

l 

l 

i 

— 

0.99 

1.02 


6 

1 

1 

1 

1 

— 

1.03 

i 

1 

l 

— 

0.97 

1.03 

761 X | 

1 

l 

1.18 

1.09 

l 

0.98 

1.10 

1.11 

GFRP | 

8 plies | 

2 

l 

l 

i 

1.20 

1.08 

l 

l 

| 

1.00 

1.02 

1.05 

10) sym | 

3 

i 

1 

i 

1.20 

1.08 

1 

1 

l 

0.98 

1.05 

1.08 

L/t = | 

4 

1 

1 

1.08 

0.99 

1 

1 

0.99 

0.94 

0.96 

109.2 | 

5 

l 

1 

1.10 

1.02 

1 

l 

I 

0.99 

0.98 

1.02 


6 

1 

l 

1 

1 

1.14 

1.03 

1 

l 

1 

0.99 

0.95 

1.01 

765 j 

1 

1 

— 

1.09 

1 

— 

0.82 

0.83 

GFRP j 

8 plies | 

2 

1 

1 

i 

— 

0.88 

l 

1 

1 

— 

0.65 

0.67 

C45/-45/ j 

3 

l 

1 

— 

0.95 

1 

l 

— 

0.68 

0.70 

4S/-45] j 
sym | 

4 

1 

l 

1 

— 

1.07 

l 

l 

1 

— 

0.78 

0.81 

L/t = | 

5 

1 

1 

— 

1.06 

1 

1 

— 

0.78 

0.81 

159.0 | 

i 1 

6 

l 

l 

1 

■ 

— 

1.02 

i 

l 

1 

l 

— 

0.72 

1 

0.77 

1 
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Table 4. (Continued) 


1 

Plate | 

Mode 

1 

l 

Natural 

Frequency Ratios, (computed / measured) 

Number | 

Number | 








and | 


l 

CLT 


Lin 

l 

Xiao | 

NASTRAN 

NASTRAN 

pa rams . | 


l 



et. al 


et. al. | 

12-by-12 

12- by- 12 

| 


l 

ref. 

12 | 

ref. 11 

l 

ref. 12 j 
1 

lumped 

consis. 

1 

769 | 

1 

l 

— 


0.95 

I 


0.93 

0.95 

GFRP | 


l 




l 

i 



8 plies | 
I 

2 

l 

| 

— 


1.08 

I 

l 

i 

1.03 

1.06 

1 

[0/90/45/| 

3 

1 

1 

— 


0.98 

i 

l 

1 

0.92 

0.94 

-45] sym | 


1 




l 

l 



i 

4 

l 

— 


1.01 

l 


0.96 

0.99 

L/t = | 


1 




I 

l 



163.6 j 
1 

5 

1 

i 

— 


1.04 

i 

1 

l 

0.99 

1.03 

1 

1 

J 

6 

1 

1 

l 

— 


1.01 

1 

i 

I 

i 

| 

0.94 

0.99 

771 1 

1 

i 

— 


1.20 

I 

| 

1.19 

1.20 

GFRP | 


1 




i 

i 



12 plies | 
j 

2 

1 

| 

— 


1.17 

I 

l 

i 

1.13 

1.16 

1 

[0/ -60/60 | 

3 

1 

1 

— 


0.98 

i 

I 

1 

0.93 

0.96 

/0/-60/60] | 


i 




i 

I 



sym | 

1 

4 

l 

1 

— 


1.06 

I 

1 

i 

1.02 

1.06 

L/t = j 

5 

1 

1 

— 


1.07 

1 

l 

1 

1.03 

1.07 

97.0 | 


l 




I 

1 



1 

l 

I-. 

6 

1 

l 



1.03 

I 

I 

I 

.. . . 1 

0.95 

1.01 

1 

762 | 

1 

1 

l 

— 


1.03 

I 


1.01 

1.02 

CFRP | 


1 




i 

i 



8 plies | 
1 

2 

l 

i 

— 


1.10 

l 


1.02 

1.05 

1 

[0] sym | 
1 

3 

1 

1 

1 

— 


1.06 

1 

I 

1 

1 

t 

0.99 

1.04 

L/t = | 

4 

I 

l 

— 


1.11 

1 

I 

1 

1.01 

1.08 

112.7 | 


1 




I 

1 



l 

i 

5 

1 

i 

— 


1.10 

l 

1 

1 

1.00 

1.07 

l 

1 

1 

- 1-- 

6 

1 

l 

1 

1--. 

— 

1 

l 

1.03 

1 

i 

i 

1 

I 

1.00 

i 

1.03 

1 1 l 1 1 1 1 
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Table 4. (Continued) 


1 

Plate | 
Member | 
and ] 

params. | 

I 

i 

Mode 

Number 

Natural 

CLT 

ref. 12 

Frequency Ratios, (computed / measured) | 

Lin | Xiao | NASTRAN | NASTRAN j 
et. al. j et. al. j 12-by-12 j 12-by-12 | 
ref. 11 j ref. 12 | lumped | consis. | 

764 | 

CFRP j 

1 

— 

1 

0.84 

1 

1 

0.81 

1 

0.82 

8 plies | 

i 

2 

— 

0.97 


0.95 

0.97 

1 

[0/90/0/ | 
90] sym | 

3 

— 

0.97 

— 

0.93 

0.96 

1 

L/t = | 

4 

- — 

0.99 


0.97 

0.99 

110.6 | 

5 

— 

1.00 

— 

0.97 

1.00 

l 

I 

I 

1 

6 

— 

0.98 

— 

0.92 

0.97 



770 | 

1 

1.14 

i.ii 

1.08 

1.10 

1.12 

1 CFRP 1 | | I 1 1 i 

| 8 plies | 

2 

1.16 

i.ii 

1.06 

1.09 

1.12 | 
1 

I i 

| [0/90/45 | 

3 

1.15 

1.09 

1.00 

1.06 

1.10 

1 /-45] sym | | I 1 1 1 1 

i i 

4 

1.07 

1.00 

1.00 

0.98 

1.01 

I L/t = | | 1 1 ! ! 1 

1 132.7 | 

5 

1.15 

1.08 

0.99 

1.06 

1.09 

I I 

I I 

I I 

i j 

6 

1.11 

1.03 

1.00 

0.99 

1.04 

772 j 

1 

1.07 

1.05 

1.00 

1.04 

1.05 

1 CFRP 1 1 | 1 1 1 

j 12 plies | 

• i 

2 

| 1.05 

| 1.03 

1.00 

| 0.99 

| 1.02 

i i 

| [0/-60/60| 

3 

| 1.07 

1.04 

1.01 

| 1.01 

1 

| 1.04 

1/0/ -60/60] I 1 1 1 1 « 

| sym | 

4 

1.10 

1 

1.06 

1 1.00 
1 

1.02 

| 1.06 
| 

tl 

-J 

5 

1.09 

1.05 

1 

1 1.00 

1.02 

1.05 

1 106.7 1 1 1 1 1 1 

1 1 
1 1 

6 

| 1.09 

1 

| 1.03 

1 

1 0.99 

1 

| 0.97 

1 

| 1.03 

1 
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Table 5. Root -mean- square errors in NASTRAN-computed natural 
frequencies for ten composite plates for uncorrected 
and corrected effective lamina elastic moduli 


Plate 

Number 


RMS error in NASTRAN -predicted frequency 



l 

1 

With uncorrected 
moduli (percent) 

| With corrected 

| moduli (percent) 

734 

1 

5.3 

1 2.9 

761 L 

1 

8.8 

1 5.5 

761 X 

1 

10.7 

1 5.7 

765 

l 

34.4 

| 26.7 

769 

l 

15.6 

1 5.4 

771 

1 

9.7 

j 10.1 

762 

1 

2.1 

| 0.6 

764 

l 

11.9 

j 9.6 

770 

1 

7.1 

j 6.6 

772 

l 

9.6 

1 2.5 


Table 6. Comparison of absolute measured and predicted 

composite plate natural frequencies of Lin et. al. 
with NASTRAN CQUAD4 computations 


Plate 


Mode 


Natural Frequencies <Hz) 


Number | 
and | 

pa rams. | 

l 

Number 

I 

| Lin et. al., 

j ref. 10-12 

| measured | computed 

i i . 

| NASTRAN 
j 12-by-12 
| limped 

NASTRAN 
12-by-12 
cons is. 

1 

734 | 

GFRP | 

1 

i l 

1 62.2 | 

66.4 

| 64.2 

65.0 

8 plies | 

I 

2 

1 «1.4 1 

131.6 

j 123.8 

127.4 

1 

[0/90/0/ | 
90] sym | 

3 

1 159.2 | 

164.5 

j 156.4 

160.9 

1 

L/t = | 

4 

| 180.5 | 

189.8 

j 177.9 

184.0 

110.7 j 
1 

5 

| 200.1 | 

208.9 

| 199.2 

206.0 

1 

1 

1 

6 

1 326.7 1 

347.2 

| 321.2 

338.4 
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Table 6. (Continued) 


1 1 

| Plate | Mode 

| Number | Number 

1 and 1 
[ pa rams . | 

l I 

■ i 

Natural Frequencies (Hz) 

Lin et. al., | NASTRAN | NASTRAN 

ref. 10-12 j 12-by-12 | 12-by-12 

measured | computed | limped | consis. 

1 1 1 


I r 

1 761 L I 



1 



78.1 

88.1 

86.9 

88.0 


j GFRP | 







| 8 plies | 

■ i 

2 

131.2 

130.7 

127.1 

130.7 


1 i 

| [0] sym I 

3 

211.5 

222.2 

215.0 

222.2 


1 1 
| L/t = | 

4 

246.0 

246.1 

232.1 

238.7 


I HI* 4 I 







1 1 

5 

287.1 

297.8 

284.1 

294.2 


1 1 

1 1 

1 1 

1 ii 

6 

362.6 

374.4 

352.9 

375.0 


| 761 X | 

1 

57.2 

62.5 

62.9 

63.7 


| GFRP | 







| 8 plies | 

■ i 

2 

90.3 

97.4 

92.5 

95.2 


1 1 

| CO) sym | 

i i 

3 

148.7 

160.5 

155.9 

161.2 


1 1 
| L/t * | 

4 

181.6 

180.2 

170.4 

175.2 


| 109.2 j 







I i 

5 

211.2 

215.9 

207.5 

214.9 


I I 

I I 

I i 

■ l 

6 

270.5 

278.8 

256.9 

272.9 


j 765 j 

1 

84.0 

91.3 

68.7 

69.6 


j GFRP j 







| 8 plies | 

i i 

2 

114.0 

99.9 

74.4 

76.6 


1 1 
| C45/-45/ | 

3 

157.0 

149.5 

107.1 

110.0 


j 45/-45J j 







| sym | 

i i 

4 

199.3 

212.6 

156.3 

161.6 


1 1 
| L/t = | 

5 

213.4 

226.7 

167.2 

173.4 


| 159.0 j 







1 1 

1 1 

• i 

6 

346.6 

353.5 

i 

249.6 

265.5 

.1 

■ -I 
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Table 6. (Continued) 


Plate 

Number 


Mode 

Number 


Natural Frequencies (Hz) 


and [ 


l 

L 

n et. al . , 

l 

NASTRAN 

l 

NASTRAN 

params. | 


1 

ref. 

10-12 

l 

12-by-12 | 

12-by-12 

1 


l 

measured | 

i 

computed | 

i 

lumped 

I 

I 

cons i s . 

1 

769 | 

1 

1 

1 

58.2 

i 

55.5 

i 

1 

54.4 

I 

55.1 

GFRP j 


1 


I 


l 


I 


8 pi ies | 
l 

2 

l 

| 

91.6 

I 

1 

99.0 

1 

i 

94.5 

I 

97.2 

I 

[0/90/45/1 

3 

1 

l 

125.5 

i 

I 

123.0 

i 

l 

115.2 

I 

I 

118.4 

-45]sym | 


1 


I 


l 


I 


l 

4 

l 

150.4 

I 

151.3 

l 

143.8 

l 

148.7 

L/t = | 


1 


i 


l 


i 


163.6 | 

1 

5 

1 

| 

156.8 

I 

1 

163.0 

l 

155.7 

I 

i 

161.1 

1 

l 

l 

i 

6 

1 

i 

l 

277.3 

1 

i 

i 

279.4 

i 

l 

l 

260.5 

i 

i 

l 

1 

274.4 

771 1 

1 

1 

l 

90.4 

I 

108.2 

1 

107.4 

l 

108.8 

GFRP j 


i 


I 


l 


i 


12 pi ies | 
| 

2 

1 

| 

144.7 

I 

1 

168.6 

l 

i 

163.6 

I 

168.4 

1 

[0/ - 60/60 | 

3 

1 

1 

222.3 

1 

i 

218.6 

i 

1 

206.6 

I 

i 

212.4 

/0/-60/60) | 


1 


i 


1 


i 


sym | 

I 

4 

i 

1 

264.1 

i 

1 

280.2 

l 

1 

270.6 

i 

i 

279.8 

L/t = | 

5 

1 

l 

281.1 

1 

l 

301.0 

i 

1 

290.9 

i 

1 

300.9 

97.0 j 


l 


l 


1 


l 


1 

I 

j .. 

6 

1 

1 

492.6 

I 

I 

505.2 

1 

1 

469.7 

l 

l 

1 _ 

499.3 

1 

762 | 

1 

' 1 " 

1 

81.5 

i 

83.6 

l 

82.2 

I 

83.3 

CFRP | 


i 


i 


1 


i 


8 pi ies | 

I 

2 

1 

i 

107.4 

I 

1 

118.4 

1 

1 

109.4 

I 

i 

112.5 

1 

[0] sym | 
1 

3 

1 

l 

| 

196.6 

1 

I 

1 

207.8 

1 

1 

1 

198.4 

i 

I 

t 

205.1 

L/t = | 

4 

1 

l 

295.5 

1 

I 

329.4 

1 

l 

299.2 

i 

1 

318.2 

1 112,7 | 


l 


l 


l 


I 


l 

1 

5 

l 

| 

382.5 

I 

j 

419.8 

l 

1 

383.6 

l 

i 

410.2 

1 

l 

l 

I... 

6 

1 

l 

1 

. 1 .. 

531.0 

1 

I 

i 

546.9 

1 

l 

l 

530.6 

i 

I 

I 

i 

545.9 

I 

' I 1 1 1 1 
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Table 6. (Continued) 


Plate 

Number 


Mode 

Number 


Natural Frequencies (Hz) 


and | 

pa rams. | 

I 

I 


Lin et. al. f 
ref. 10-12 
measured | computed 

i 

NASTRAN 

12-by-12 

lumped 

NASTRAN 

12-by-12 

consis. 

764 | 

1 

i 

68.9 | 

l 

58.1 

55.5 

56.2 

CFRP 1 | 1 1 1 1 

8 plies | 

2 

218.9 

213.3 

206.9 

212.8 

i 

[0/90/0/ | 

3 

251.2 

243.5 

233.3 

241.5 

90] sym | 

i 

4 

305.4 

302.5 

| 294.9 

303.4 

L/t = | | I 1 1 

110.6 | 

5 

323.5 

324.2 

312.7 

323.7 

l 

l 

l 

6 

452.5 

441.6 

414.5 

437.8 

770 j 

1 

77.8 

86.3 

85.9 

87.0 

CFRP | | 1 i 1 

8 plies | 

2 

| 202.7 

| 224.5 

| 220.6 

| 226.9 

i 

(0/90/45 | 

3 

| 258.0 

j 280.4 

| 273.9 

j 283.4 

/-45] sym | 

I 

4 

| 298.7 

| 298.8 

| 293.6 

j 302.0 

1 L/t = 1 | 1 1 1 

132.7 | 

5 

| 322.0 

| 348.4 

| 340.6 

| 352.4 

l 

l 

l 

6 

| 496.7 

j 512.2 

| 490.0 

| 517.3 

1" 

772 | 

1 

| 156.6 

| 165.2 

| 162.7 

j 165.0 

1 CFRP 1 | I 1 1 

| 12 plies | 

2 

| 272.0 

i 

| 279.1 

i 

| 269.9 

I 

| 277.6 

1 

l 1 

| [0/-60/60| 

3 

1 

| 372.3 

1 

| 387.8 

1 

| 376.9 

1 

| 387.6 

|/0/- 60/60] | 
| sym | 

4 

I 

| 407.8 

1 

l 

| 432.6 

1 

l 

| 416.7 

1 

1 

| 431.1 

1 

| 1 
L/t = | 

5 

1 

| 486.1 

1 

| 511.4 

1 

| 494.7 

1 

| 512.0 

1 106.7 1 | | 1 1 

1 l 

1 l 

6 

| 779.0 

.1 

| 800.4 

.1 

| 752.6 

.1 

| 799.4 

.1 
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Table 7. Mode-by-mode tabulation of NASTRAN CQUAD4 idealization 
frequency prediction errors, with averaged percentage 
errors and eigenmode descriptions 



Plate 

l 

l 

Mode 

i 

j Percentage errors 

l 

| E i genmode 



N uribe r 

I 

Number | 



| description 



and 

l 


| Lin et. 

NASTRAN 

| NASTRAN 

I 



params. 

l 


j al. 

12-by-12 

| 12-by-12 | 




I 


i 

-- 1 

lumped 

consis. 

i 



734 

I 

1 

1 

1 + 7 

+ 3 

+ 5 

| shear and flexure 



GFRP 

l 


I 



i 



8 plies | 

2 

1 0 

- 6 

- 3 

| 2-noded beam abt. 




i 


i 



| weak axis 



[0/90/0/ | 

3 

1 + 3 

- 2 

+ 1 

| 2-noded beam abt. 



90] sym 



i 



j strong axis 





4 

1 + 5 

- 1 

+ 3 

| 2-noded beam abt. 



L/t = 



I 



| weak axis + shear 



110.7 


5 

1 + 4 

0 

+ 3 

| 2-noded beam abt. 






I 



| strong axis + shear 1 




6 

1 + 6 

* 2 

+ 4 

fundamental 






i 



plate flexure 





avg. 

j ♦ 4.2 

- 1.3 

+ 2.2 




761 L 


1 

j ♦ 13 

+ 11 

+ 13 

shear and flexure 



GFRP 



I 






8 plies 


2 

1 0 

- 3 

0 

2-noded beam abt. 






l 



weak axis 



[0] sym 


3 

1 + 5 

+ 2 

+ 5 

2-noded beam abt. 






i 



weak axis + shear 



L/t = 


4 

1 0 

- 6 

- 3 

2-noded beam abt. 



111.4 



i 



strong axis 





5 

1 + * 

- 1 

♦ 2 

2-noded beam abt. 






1 



strong axis + shear 





6 

1 + 3 

- 3 

+ 3 

3-noded beam abt. 






| 



weak axis 





avg. 

j + 4.2 

0 

+ 3.3 




761X 


1 

1 + 9 

+ 10 

+ 11 

shear and flexure 



GFRP 



i 






8 plies 


2 

I + 8 

+ 2 

♦ 5 

2-noded beam abt. 






l 



weak axis 



[0] sym 


3 

l + 8 

+ 5 

+ 8 

2-noded beam abt. 






i 



weak axis + shear 



L/t = 


4 

I ■ 1 

- 6 

- 4 

2-noded beam abt. 



109.2 



I 



strong axis 





5 

1 + 2 

- 2 

+ 2 

2-noded beam abt. 






1 



strong axis + shear 





6 

| + 3 

- 5 

♦ 1 

3-noded beam abt. 






. 1 j 



weak axis 





avg. 

1 1 
1 ♦ 1 

'1 1 

+ 0.7 

+ 3.8 
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Table 7. (Continued) 


1 1 
| Plate | 

Mode 

I 

Percentage errors ] 

E i genmode 


| Number | 

Number 




description 


1 and 1 


Lin et. 

NASTRAN 

NASTRAN 



| pa rams . | 


al. 

1Z-by-12 

12-by-12 



1 1 

l t 



limped 

consis. 



1 765 1 

1 

+ 9 

- 18 

- 17 

plate flexural 


j GFRP | 







| 8 plies | 

i i 

2 

- 12 

- 35 

- 33 

plate flexural 


1 1 
j C45/-45/ | 

3 

- 5 

- 32 

- 30 

plate flexural 


j 45/-451 | 







| sym | 

i i 

4 

♦ 7 

- 22 

- 19 

plate flexural 


1 1 
| L/t = | 

5 

♦ 6 

- 22 

- 19 

plate flexural 


| 159.0 | 







1 1 

l l 

l 1 

6 

♦ 2 

- 28 

- 23 

plate flexural 


l l 

avg. 

+ 1.2 

- 26.2 

- 23.5 



1 769 1 

1 

- 5 

- 7 

- 5 

shear and flexure 


j GFRP | 







| 8 plies | 

2 

♦ 8 

+ 3 

♦ 6 

2-noded beam abt. 


I I 





weak axis 


j [0/90/45/ j 

3 

- 2 

- 8 

- 6 

2-noded beam abt. 


| -453 sym | 





strong axis 


1 1 

4 

♦ 1 

- 4 

- 1 

plate flexural 


j t/t = | 







j 163.6 j 

5 

♦ 4 

- 1 

+ 3 

plate flexural 


I i 
i i 
I i 

i 1 ■ 

6 

♦ 1 

- 6 

- 1 

plate flexural 


i i 

i i 

avg. 

♦ 1.2 

- 3.8 

- 0.7 



1 771 1 

1 

+ 20 

+ 19 

+ 20 

shear and flexure 


| GFRP | 







| 12 plies | 

2 

+ 17 

+ 13 

+ 16 

2-noded beam abt. 


| | 





weak axis 


| [0/ - 60/60 | 

3 

- 2 

- 7 

- 4 

2-noded beam abt. 


j/0/-60/603 j 





strong axis 


| sym | 

■ i 

4 

♦ 6 

+ 2 

+ 6 

plate flexural 


1 1 
| L/t = | 

5 

+ 7 

+ 3 

♦ 7 

plate flexural 


I 97 -° 1 







| | 

6 

♦ 3 

- 5 

♦ 1 

3-noded beam abt. 


1 1 





weak axis 


1 i 

l 1 

i i 

avg. 

♦ 8.5 

i 

+ 4.2 

+ 7.7 

i 


-- 1 
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Table 7. (Continued) 


1 

Plate | 

Mode 

I 

| Percentage errors 

E i genmode 

Nunber | 

Number 




description 

and | 


Lin et. 

| NASTRAN 

| NASTRAN 


pa rams. | 


al . 

| 12-by-12 | 12-by-12 


l 



lumped 

| consis. 


762 j 

1 

+ 3 

♦ 1 

j + 2 

shear and flexure 

CFRP | 




l 


8 pi ies | 

2 

+ 10 

+ 2 

1 + 5 

2-noded beam abt . 

l 




l 

weak axis 

[0] sym | 

3 

+ 6 

- 1 

1 + 4 

2-noded beam abt. 

I 




l 

weak axis + shear 

L/t = | 

4 

+ 11 

+ 1 

1 + 8 

3-noded beam abt. 

112.7 | 




l 

weak axis 

1 

5 

+ 10 

0 

1 * 7 

3-noded beam abt. 

l 




I 

weak axis + shear 

1 

6 

+ 3 

0 

1 + 3 

2-noded beam abt. 

1 




| 

strong axis 

| 

avg. 

+ 7.2 

+ 0.5 

j + 4.8 


1 

764 | 

1 

- 16 

- 19 

j - 18 

shear and flexure 

CFRP | 




I 


8 plies | 

2 

- 3 

- 5 

1 ' 3 

2-noded beam abt. 

1 




i 

weak axis 

[0/90/0/ | 

3 

* 3 

- 7 

1 - 4 

2-noded beam abt. 

90] sym | 




l 

weak axis + shear 

1 

4 

- 1 

- 3 

1 - 1 

2-noded beam abt. 

L/t = | 




I 

strong axis 

110.6 j 

5 

0 

- 3 

1 0 

2-noded beam abt. 

1 




l 

strong axis + shear 

1 

6 

- 2 

- 8 

l - 3 

fundamental 

1 




| 

plate flexure 

1 ‘ 
1 

avg. 

- 4.2 

- 7.5 

j - 4.8 


1 

770 | 

1 

+ 11 

+ 10 

j + 12 

shear and flexure 

CFRP | 




1 


8 pi ies | 

2 

+ 11 

♦ 9 

1 + 12 

2-noded beam abt. 

l 




1 

weak axis 

[0/90/45 | 

3 

1 + 9 

♦ 6 

| + 10 

plate flexural 

/-45] sym | 




1 


1 

4 

0 

- 2 

1 + 1 

2-noded beam abt. 

L/t = | 




1 

strong axis 

132.7 j 
1 

5 

+ 8 

♦ 6 

1 + 9 
| 

plate flexural 

1 

1 

1 

6 

+ 3 

- 1 

1 

1 + 4 
1 

plate flexural 

1" 

1 

I- 

avg. 

+ 7.0 

+ 4.7 

j+8.0 | 

1 ! 
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Table 7. (Continued) 


Plate | 
Nunber | 
and | 

pa rams. | 

Mode 

Number 

Percentage errors 

Lin et. | NASTRAN | NASTRAN 

al. | 12-by-12 | 12-by-12 

| lumped | consis. 

772 | 

1 

♦ 5 

1 

1 + 4 

l 

1 4 5 

CFRP j 

12 plies | 

2 

♦ 3 

1 

1 ‘ 1 
1 

| + 2 

t0/-60/60| 

3 

♦ 4 

1 

1 + 1 

1 + 4 

/0/-60/60] j 

sym | 

4 

+ 6 

1 

1 4 2 
1 

1 * 6 

L/t = j 

5 

♦ 5 

1 + 2 

1 + 5 

106.7 j 

6 

+ 3 

l 

I ■ 3 
| 

1 + 3 


avg. 

♦ 4.3 

| + 0.8 

j + 4.2 


Eigenmode 

description 


shear and flexure 

2-noded beam abt. 

weak axis 
2-noded beam abt. 
strong axis 
plate flexural 

plate flexural 

3*noded beam abt. 
weak axis 
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Notes: Positive X-axis: fiber orientation angle of 0 degrees 
Positive Y-axis: fiber orientation angle of 90 degrees 


8-noded rectangles 4-noded CQUAD4's 



(i) Lin and Xiao et. al. (ii) NASTRAN CQUAD4 

idealization idealization 


Figure 2. Plate element meshes used by Lin and 
Xiao et. al. and in NASTRAN analyses 
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FLHTE M3. 73* OF XIRQ,LIN. JU C 0 / 90 / 0/30 3 S GFRP < I E-BV- 1 2 ) 
SOLUTION 3 - EICB^MRLUE flNHLVSIS 
CSC 6;W B;EICU £.£SBBE+OS3 


Fig 3b. Plate 734, Mode 2 

Figure 3. NASTRAN-computed eigenmodes for plate 734 
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CSC C.*MS 9;EIG*/ 1.0C51E+OC3 


Fig. 3C. Plate 734, Mode 3 



CSC 10:E!GU 1 .40S4E+0C3 


Fig. 3d. Plate 734, Mode 4 


Figure 3. (Continued) 
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Fig. 3e. Plate 734, Mode 5 



Fig. 3f . Plate 734, Mode 6 

Figure 3 . (Continued) 
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PLffre T* VF L!N.NI,FWB (ISM) C <0>83 gfvf ue-^-i5> 
SOLUTION a - EICB^HLUE flNRLVSIS 
CSC T;rr« t;EICL' e.S643E*fOS3 

Fig. 4a . Plate 761L, Mode 1 



Fig 4b. Plate 761L, Mode 2 


Figure 4 . 


NASTRAN-computed eigenmodes for plate 761L 




Fig. 4f. Plate 761L, Mode 6 

Figure 4. (Continued) 
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Plate 765, Mode 2 


SOLUTION 3 - EIGETM=CUE RWLVSIS 
CSC B;f-fS B ; E I CU 1.G633E+OS3 


Fig 5b. 


Figure 5. 


NASTRAN-computed eigenmodes for plate 765 





SOLUTION a - E I OE7*/FUUE RTflL^IS 
CSC B;i-H B.’EICU 1.S4TTE+OS3 


Fig. 5c. Plate 765 # Mode 3 



CSC IO.‘Ki lOJElQ^ 7.7B70E+0S3 


Fig. 5d. Plate 765, Mode 4 


Figure 5. (Continued) 
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Fig. 5f . Plate 765, Mode 6 

Figure 5. (Continued) 
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Fig. 6d. Plate 769, Mode 4 


Figure 6. (Continued) 
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PLHTE TG9 GF Lm.NMOtC < I9B4> C0.90.«,-«3S GFFR ( 1 E-&t'-I 2> 
sotxrricM 3 - eigenumjle rtf*. vs is 
CSC u:m 1 1 :EI»/ T.4TME+0S3 


Fig. 6e. Plate 769, Mode 5 



Fig. 6f. 


Plate 769, Mode 6 


Figure 6. (Continued) 
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esc e.-m e.Eicu b.thi ie-kjs: 

Fig 7b. p; 

Figure 7 . NASTRAN-comp' 
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Figure 7. (Continued) 
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Fig. 7f. Plate 771, Mode 6 

Figure 7. (Continued) 
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Fig. 8e. Plate 762, Node 5 



Fig. 8f . Plate 762, Mode 6 


Figure 8. (Continued) 
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N92-24330 


A CASE OF POOR SUBSTRUCTURE DIAGNOSTICS 
Thomas G. Butler 
BUTLER ANALYSES 


Substructuring is a powerful tool. As with any powerful 
tool the options for managing a job are legion. On the other 
hand, the NASTRAN Manuals in the Substructuring area are all 
geared toward instant success, but the solution paths are fraught 
with many traps for human error. Thus, the probability of suf- 
fering a fatal abort is high. In such circumstances, the neces- 
sity for diagnostics that are user friendly is paramount. This 
paper is written in the spirit of improving the diagnostics as 
well as the documentation in one area where the author felt he 
was backed into a blind corner as a result of his having com- 
mitted a data oversight. This topic will be aired by referring 
to an analysis of a particular structure. 

The structure, under discussion, used a number of local 
coordinate systems that simplified the preparation of input data. 
The principal features of this problem are introduced by refer- 
ence to a series of figures. 

Figure 1 illustrates a PILOT model of the basic component 
substructure of a full scale structure. This pilot model 
was used to explore the error that developed in the true 
structure. In preparation for the investigation into the 
difficulty that was encountered during a "COMBINE" operation, 
the pilot basic was cloned 4 times into CLONA, CLONB, CLONC 
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and CLOND. 

Figure 2 tabulates the bulk data lor the 5 coordinate systems 
that were used in the basic component. Coordinate system "5" 
is cylindrical and was used for its core. Coordinate systems 
”50, 60, 70 & 80" are rectangular and were used for the four 
arms with their local X axes pointing outward at zero de- 
grees, 90, 180 and 270 respectively. Each clone retained its 
own copy of the set of five local coordinate systems. Thus, 
the Substructure Operating File (SOF) at this point had a 
complement of 5x5= 25 coordinate systems to catalog. The 
multiplicity of coordinate systems was at the root of the 
fatal error which erupted. 

Figure 3 illustrates two separate "COMBINE" operations 
amongst the substructures. In the first "COMBINE", point 51 
of P/S CLONC joins with point 71 of P/S CLOND. In the second 
"COMBINE", point 61 of P/S 1 BASE joins with point 81 of P/S 
CLONA, while point 61 of P/S CLONA joins with point 81 of P/S 
CLONB. 

During the subsequent linking of substructures , the points that 
were combined each had their own local coordinate systems. Well 
this doesn't seem to be a problem, because NASTRAN has a wonder- 
ful module called CSTM (Coordinate System Transformation Matrix) 
which keeps track of all transformations amongst a host of coor- 
dinate systems. So the user is disarmed into thinking that 


1. The abbreviation P/S, meaning pseudo-structure, is used as a 
generic term for any number of different kinds of substructures, 
basic, or clones, or condensations, or combined. 
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NASTRAN can handle anything dealing with coordinates. This was 
especially true in this case, because, just prior to the abort 
being described here, a mistake in translating one of the cloned 
structures was corrected in response to a diagnostic message that 
declared that points, which were targeted to be joined, were not 
within tolerance. The error was that one of the translations, 
defined on a "TRANS" card was off by an eighth of an inch. After 
the correction a message was issued declaring that all points in 
the "COMBINE" operation were within tolerance. So the reaction 
to a subsequent message to the effect that the local coordinate 
systems were incompatible seemed ridiculous, because NASTRAN had 
no difficulty in locating the points in space and in pronouncing 
that they were within tolerance with the coordinate systems that 
were corrected. 

As it turned out there are a number of different coor- 
dinate systmes that have to be dealt with here, and the "TRANS" 
set that was just corrected - though at first suspected - was not 
at the nub of the problem. The problem arises not in the align- 
ment, which the TRANS coordinates deal with, but in the sub- 
sequent mating, which depends of the local DISPLACEMENT coor- 
dinates of points that are being brought together. 

As a matter of general substructuring principle, when a 
group of substructures is assembled, any place where parts are 
linked can involve contributions from 2 or more individuals. At 
any such place the set of points are merged into a resultant 
single point. What is not told in the manuals is that the resul- 
tant point needs to refer to just one coordinate system. If all 
of the merging points refer to a common coordinate system, there 
is no problem. But, when each point has its own local displace- 
ment coordinate system, NASTRAN aborts and issues a message #6528 
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saying that incompatible local coordinate systems have been 
found. But if the user thinks that the problem has just been 
corrected, the characterization of the coordinates in message 
#6528 as being “incompatible” doesn't make sense and he becomes 
convinced that there must be a bug in NASTRAN and the user is to 
be absolved of blame. His certitude of blamelessness is further 
reenforced by the details that are supplied with the diagnostic 
message. The text of the complete message, shown in Figure 4, 
refers to local coordinate systems 1 and 10. But if you look at 
Figure 2, you can verify that no coordinate system was numbered 1 
or 10. This seems to further corroborate that NASTRAN got some 
tables mixed up and is in need of having a bug straightened out. 

Gordon Chan of the UNISYS Support Group came to my rescue 
and published the transformation matrices for the coordinate 
systems that were involved. The reason that NASTRAN aborted was 
uut because ui. an error in the code. It deliberately compared 
the local coordinates at the combining point and found that one 
pair of signs was aligned while the signs of two other axes were 
of opposite in sign. The message referring to coordinate 
systems resulted from a partially completed execution of the 
COMBINE command. It had reassinged coordinate system ID's in 
terms of its own internal bookkeeping system, but it phrased the 
diagnostic in terms of its own scheme of ID's. Unf ormtunately , 
that part of its completed operation was never output, because of 
the abortion, so the diagnostic which was trying to be helpful 
was confusing the situation even further. However, NASTRAN 
appeared to be operating properly. 

Double checking of coordinate systems 50 through 80 found 
them to be error free. As a further check, the manual method was 
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compared to the automatic method of combining. The same diagnos- 
tic regarding incompatible local coordinates showed up in this 
automatic trial as well, but this time referred to pairs of 
coordinates with other sets of strange identifications; i.e. 2 & 
9 and 7 & 14. That diagnostic is shown in Figure 5. Finally a 
vague, misty fragment seemed to kindle in the back of my brain 
that had something to do with the data card called "GTRAN". I 
pored over the Substructure Section of Chapter 1 of the User's 
Manual to uncover a hint on the use of GTRAN. No help. Nor was 
the Theoretical Manual any assistance. Figure 1 shows that in 
the example of the manual COMBINE, points 51 and 71 refer to 
coordinate systems 50 & 70 respectively. NASTRAN finds that 
these two systems do not align with each other and so both cannot 
be allowed to represent that point after a merge. The situation 
must be reconciled and NASTRAN needs guidance from the user. The 
avenue by which the user exerts his preference is through the use 
of GTRAN. The bulk data explanation of GTRAN left many unanswer- 
ed questions. The only thing left to do was to resort to the old 
"black box" method of finding out how it behaved. GTRAN was 
tried out under its options. One option is to refer all connect- 
ing points to the overall basic system, and the other is to refer 
them to the system defined by the TRANS entry. Both worked! 
Figure 6 shows an excerpt of the output from a successful manual 
run using GTRAN. It repeats the message about points within 
tolerance, then gives the tabulation of the resulting points 
after the COMBINE operation, showing those degrees of freedom 
that were merged into a single point. This connectivity summary 
does not, however, refer to any coordinate system. Coordinate ID 
information is published subsequently in the BGSS. In this case 
the BGSS shows that it was arbitrated by referring both points to 
the "0" system (the overall basic). 


104 



SUBSTRUCTURE DIAGNOSTICS 


Figure 7 is an excerpt of a summary of connectivities for 
the automatic COMBINE case after a proper use of GTRAN. It shows 
a similar set of connections as in the manual case but amongst 
BASE, CLONA, and CLONB. 

There were many unhappy features relating to documenta- 
tion in this encounter: (1) the diagnostic itself, (2) the 
explanation of the diagnostic in Chapter 6, (3) the guide to 
modeling in Chapter 1, (4) the explanatory notes in the bulk 
data, and (5) the Theoretical Manual. It is incumbent upon the 
manuals to acquaint the user with what its needs are so that he 
can supply necessary data. But in this instance the documenta- 
tion gave NO hint of how NASTRAN operated internally, so the user 
was set adrift by a diagnostic that impugns his data as INCOMPAT- 
IBLE. For all he knew NASTRAN had some sort of internal default 
to meet the arbitration needs. Without the help of documenta- 
tion, the user must look into the code to find out what NASTRAN 
is doing in subroutine "COMBI". He does not know from the above 
documentary sources whether NASTRAN takes a default when not 
supplied with specific direction or aborts. The situation is 
this. NASTRAN first determines that the points that it is di- 
rected to link are collocated. This can be done by temporarily 
transforming all locations to the overall basic system. But now 
when it wants to trim all connecting points to a single point, it 
must assign some coordinate system to that resulting point. But 
which one? Dave Her ting and the savants that helped him with the 
architecture of SUBSTRUCTURING were aware of the problem and 
provided for it with the GTRAN card. But as is often the case 
with programming, the documentation did not coach the user into 
anticipating the need to guide NASTRAN in the assignment of a 
coordinate system to a common point. 
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Rather than overcome the obstacle with the provision of a 
GTRAN card and then to continue with the analysis of the struc- 
ture only, I chose to share this experience at the Colloquium and 
to volunteer a supplement to the documentation so that any 
subsequent user can be well quided when he encounters message 
#6528. Figure 8 shows the recommended diagnostic message. 
Figure 9 shows the recommended supplement to the "COMBINE" 
section of Chapter 1 on modeling with substructures, and in 
Chapter 6 on explanation of diagnostics. No suggestions are 
offered for the Theoretical Manual, because it is currently 
awaiting a major revision. 

I extend my deep appreciation to Gordon Chan for his help 
in unearthing this problem and for his modification of the diag- 
nostic message in the code. The new release will have the re- 
vised diagnostic message. In addition Gordon Chan added a print- 
out of the transformation matrices of the coordinate systems that 
are indicted. 

My hope is that this small effort will save future users 
much time and frustration when faced with an unsuccessful COMBINE 
operation in their substructuring work. 
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C0RD2C 

5 

0 

0.0 

0.0 

0.0 

0.0 

0.0 

10.0 

+CYLN 

+CYLN 

10.0 

0.0 

10.0 







C0RD2R 

50 

0 

10.0 

0.0 

10.0 

10.0 

-10.0 

10.0 

+RAY0 

+RAYO 

20.0 

-10.0 

10.0 







C0RD2R 

60 

0 

0.0 

10.0 

10.0 

10.0 

10.0 

10.0 

+RAY90 

+RAY90 

10.0 

20.0 

10.0 







C0RD2R 

70 

0 

-10.0 

0.0 

10.0 

-10.0 

10.0 

10.0 

+RAY18 

+RAY180 

-20.0 

10.0 

10.0 







C0RD2R 

80 

0 

0.0 

-10.0 

10.0 

-10.0 

-10.0 

10.0 

+RAY27 

+RAY270 

-10.0 

-20.0 

10.0 








Figure 2. Coordinate Systems in Component BASE 
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> > 


TRANS TRANS 

40 80 


Figure 3. Connection Diagram of Two COMBINE Operations 
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USER INFORMATION MESSAGE 6516, 

ALL MANUAL CONNECTIONS SPECIFIED ARE 
ALLOWABLE WITH RESPECT TO TOLER. 


USER FATAL MESSAGE 6528, 

INCOMPATABLE LOCAL COORDINATE SYSTEMS 
HAVE BEEN FOUND. CONNECTION OF POINTS 
IS IMPOSSIBLE, SUMMARY FOLLOWS. 


************************************** 


THE FOLLOWING MISMATCHED LOCAL COORDINATE 
SYSTEMS (CSTM) HAVE BEEN FOUND FOR 


LOCAL COORDINATE SYSTEM ID NO. 1 

PSEUDOSTRUCTURE ID NO. 1 

INTERNAL POINT NO. 2 


*************************************** 


LOCAL COORDINATE SYSTEM ID NO. 10 

PSEUDOSTRUCTURE ID NO. 2 

INTERNAL POINT NO. 14 

USER FATAL MESSAGE 6537, MODULE COMBI 
TERMINATING DUE TO ABOVE ERRORS. 


Figure 4. Diagnostic From Abort of Manual COMBINE 
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SUMMARY OF AUTOMATICALLY GENERATED CONNECTIONS 


CONNECTED 
DOF 
123456 
123456 
123456 
123456 
123456 
123456 


CONNECTION PSEUDO STRUCTURE 


CODE 

12 

12 

12 

23 

23 

23 


BASE 
7 
5 
3 
0 
0 
0 


CLONA 
15 
13 
11 
7 
5 
3 


USER FATAL MESSAGE 6528, 

INCOMPATABLE LOCAL COORDINATE SYSTEMS 
HAVE BEEN FOUND. CONNECTION OF POINTS 
IS IMPOSSIBLE, SUMMARY FOLLOWS. 


NAMES 

CLONB 

0 

0 

0 

15 

13 

11 


*************************************************** 


THE FOLLOWING MISMATCHED LOCAL COORDINATE 
SYSTEMS (CSTM) HAVE BEEN FOUND FOR 

LOCAL COORDINATE SYSTEM ID NO. 2 

PSEUDOSTRUCTURE ID NO. 1 

INTERNAL POINT NO. 5 

************************************************ 


LOCAL COORDINATE SYSTEM ID NO. 9 

PSEUDOSTRUCTURE ID NO. 2 
INTERNAL POINT NO. 13 

************************************************* 


THE FOLLOWING MISMATCHED LOCAL COORDINATE 
SYSTEMS (CSTM) HAVE BEEN FOUND FOR 

LOCAL COORDINATE SYSTEM ID NO. 7 

PSEUDOSTRUCTURE ID NO. 2 
INTERNAL POINT NO. 5 

************************************************ 

LOCAL COORDINATE SYSTEM ID NO. 14 

PSEUDOSTRUCTURE ID NO. 3 
INTERNAL POINT NO. 13 

USER FATAL MESSAGE 6537, MODULE COMBI 
TERMINATING DUE TO ABOVE ERRORS. 


Figure 5. Diagnostic From Abort of Automatic COMBINE 
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USER INFORMATION MESSAGE 6516, 

ALL MANUAL CONNECTIONS SPECIFIED ARE 
ALLOWABLE WITH RESPECT TO TOLER. 


SUMMARY OF PSEUDOSTRUCTURE CONNECTIVITIES 


INTERNAL 
POINT NO 

INTERNAL 
DOF NO 

DEGREES OF 
FREEDOM 

PSEUDOSTRUCTURE NAMES 
CLONC CLOND 

16 

89 

123456 

CLONC 

72 


17 

95 

13 

CLONC 

51 

CLOND 

71 

18 

97 

123456 


CLOND 

51 


Figure 6 

Summary of Connectivities After GTRAN Use in Manual COMBINE 
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SUMMARY OF AUTOMATICALLY GENERATED CONNECTIONS 


CONNECTED 

DOF 

123456 

123456 

123456 

123456 

123456 

123456 


CONNECTION 

CODE 

12 

12 

12 

23 

23 

23 


PSEUDOSTRUCTURE NAMES 


BASE CLONA CLONB 

7 15 0 

5 13 0 

3 11 0 

0 7 15 

0 5 13 

0 3 U 


SUMMARY OF PSEUDOSTRUCTURE CONNECTIVITIES 
INTERNAL INTERNAL DEGREES OF * PSEUDOSTRUCTURE 
POINT NO DOF NO FREEDOM BASE CLONA 


NAMES 

CLONB 


16 

89 

123456 

BASE 

72 



17 

95 

13 

BASfc 

61 

CLONA 

81 


18 

97 

123456 


CLONA 

52 


33 

183 

123456 


CLONA 

2 


34 

189 

13 


CLONA 

61 

CLONB 

81 

35 

191 

123456 



CLONB 

52 


Figure 7 

Summary of Connectivities After GTRAN Use in Automatic COMBINE 
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USER FATAL MESSAGE 6528 

INCOMPATABLE LOCAL COORDINATE SYSTEMS HAVE BEEN FOUND. COMPLETTION OF 
CONNECTION IS IMPOSSIBLE. SUGGEST USE OF "GTRAN". SUMMARY IN TERMS OF 
JUST-FORMED INTERNAL FREEDOMS AND INTERNAL COORDINATE SYSTEM ID'S PER THE 
EQSS & BGSS FOLLOW: 


Figure 8 


Revised Fatal Diagnostic Message 6528 



SUBSTRUCTURE DIAGNOSTICS 


USER'S MANUAL CHAPTER 1. 


ADD THE FOLLOWING TEXT TO SUPPLEMENT THE TOPIC OF ™E "COMBINE" OPERATION ON 
SUBSTRUCTURING IN THE NASTRAN USER'S MANUAL, PAGE 1.10-39 (14 LINES UP FROM 
THE BOTTOM OF THE PAGE. 


WHEN POINTS ARE ALIGNED FOR COMBINING AFTER A TRANSLATION ^D/OR ROTATION OF 
COMPONENTS THEY BECOME A SINGLE POINT UPON LINKING. IF THE POINTS ABOUT TO 
BE CONNECTED REFER TO DIFFERENT LOCAL COORDINATE SYSTEMS, THE SUBSTRUCTTOE 
ROUTINE "COMBI" DOES NOT IMPOSE A DEFAULT CORRDINATE SYSTEM FOR THE POINT. 
SUCH^A SITUATION MUST BE ANTICIPATED BY THE ANALYST TO AVOID A FATAL ABORTION. 
M CAN ASSIGN A DISPLACEMENT COORDI ANTE SYSTEM TO THB RESULTINC POINT 

THROUGH THE USE OF THE GTRAN CARD. IT OFFERS 3 OPTIONS: ( 1 2 J^™2 R ™I2rJ HE 

OVERALL BASIC SYSTEM, (2) NO TRANSFORMATION, AND (3) TRANSFORM TO THE COORDI 
NATE SYSTEM WHICH WAS DEFINED ON THE SELECTED "TRANS" CARD. 


USER'S MANUAL CHAPTER 6. 

ADD THE FOLLOWING TEXT AFTER THE FIRST SENTENCE OF DIAGNOSTIC MESSAGE 6528. 


EACH POINT IS CARRYING ITS OWN LOCAL COORDINATE SYSTEM INTO THE "COMBINE' D 
POINT AND THEY HAVE BEEN FOUND TO BE DIFFERENTLY ALIGNED; I.E. 

INCOMP AT AB LE . THE USER IS REQUIRED TO ARBITRATE BETWEEN THE COMPETING 
LOCAL COORDINATE SYSTEMS. HE IS ADVISED TO CONSIDER USING ONE c°5.mTufS AL 
"GTRAN" CARDS. (SEE PAGE 1.10-39 OF THE USER'S MANUAL.) HE IS FURTHER 
ADVISED TO "DESTROY" THE PSEUDO- STRUCTURE DEFINED IN THE COMBINE OPERATION 
IN ORDER TO REMOVE ANY PARTIALLY COMPLETED "COMBINE" DATA FROM THE SOF 
(SUBSTRUCTURE OPERATING FILE), BEFORE RERUNNING THE "COMBINE" OPERATION. 


Figure 9. Supplements to Documents in USER'S Manual 
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ALTERNATIVE METHODS TO MODEL FRICTIONAL 
CONTACT SURFACES USING NASTRAN 

Joseph Hoang 

GE Goverment Services, Houston, Texas 
SUMMARY 

Elongated (slotted) holes have been used extensively for the integration 
of equipment into Spacelab racks. In the past, this type of interface has been 
modelled assuming that (i) there is no slippage between contact surfaces, or 
(it) there is no load transfer in the direction of the slot. Since the contact sur- 
faces are bolted together, the contact friction provides a load path determined 
by the normal applied force (bolt preload) and the coefficient of friction. This 
paper examines three alternate methods that utilize spring elements, exter- 
nally applied couples, and stress dependent elements to model the contacted 
surfaces. Results of these methods are compared with results obtained from 
methods that use GAP elements and rigid elements. 

INTRODUCTION 

Elongated holes have been used in the design of Spacelab Experiment Equipment 
mounting provisions. This type of joint is employed where large tolerances are allowed 
in one direction of the hole for the ease of integration. A simple way to model these 
interfaces is to use RIGID elements with the assumption that two connecting grid points 
are not moving against each other when loads are applied. Another common method is to 
use RIGID elements with the degrees of freedom associated with the longitudinal direction 
of the elongated hole released. When using this method, it is assumed that the joint did 
not carry load in its slotted direction. 

Due to the assumption involved, neither method yields realistic results. GAP elements 
have been used to achieve better results. However, when large numbers of GAP elements 
are used in a complex model, an unreasonably large amount of computer processing time 
is required to solve the system and, hence, is not economical or practical. 

In this paper, three alternative methods are investigated. The first method employs 
a spring element with a spring rate equal to the maximum load at the connection divided 
by the gap length. By using a spring element at the connection, the nonlinear frictional 
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force developed at the gap is replaced by a spring force that increases proportionally with 
gap distance. The second method used an externally applied couple to represent the 
frictional forces at the two contacted surfaces. This method creates a local realization of 
frictional forces when two grid points, representing two contacted surfaces, move against 
each other. In the third method, two contacted surfaces are connected by elements with 
stress-dependent material properties. The piecewise linear static analysis rigid format is 
ultilized to solve element forces at these elements. 

A simple stowage container and four supporting columns were developed to represent 
a typical installation in a Spacelab rack. The container was integrated on supporting 
columns using the above methods. The models were run on a SUN workstation using 
CSA/NASTRAN. The results obtained from each method were then compared. 

THEORY AND BACKGROUND 

Figure la shows a typical elongated hole and Figure lb shows force vs. relative 
displacement of an elongated hole. The force vs. displacement curves of other elements 
are shown in 2 through 7. 

From inspection, the stress dependent material element is the more appropriate ele- 
ment to simulate elongated hole behavior because the GAP element, spring element and 
coupled element would generate some error. A rigid element with all three translational 
degrees of freedom coupled (R123) should be used only when the frictional force is greater 
than the force at the connection. A rigid element with degrees of freedom associated with 
the longitudinal direction of the hole released (R23) should be used when, the frictional 
force and the displacement in elongated direction are small. 


117 






1002 1004 



Fig. 8 Plot of test subject a simple box and four supported columns 
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ANALYSIS AND RESULTS 


A simple box is constructed using QUAD2 elements. It is supported by four columns 
which were modelled using BAR elements. A plot of the structure is shown in Figure 8. 
In the rear, the box is mounted on two rear columns at grid points 13, 24, 37, and 48 
using rigid elements with all three translational degrees of freedom coupled. In front, the 
box is mounted to the front left column at grid points 17 and 41 and to the front right 
column at grid points 20, and 44, using a different integrating method. Gravitational loads 
varying from 0.2g to 4.0g are applied to the model in the X-direction, which is pa.ra.11p! to 
the slotted direction. 

The material and properties of QUAD2 and BAR elements are chosen such that the 
box is heavier and more flexible than the columns. Since the columns are much stiffer than 
the box, displacement of the grid point on the column is small as compared to displacement 
of the grid point on the box. Displacement of the rear panel is closed to the displacement 
of rear columns because the rear panel interface points are mounted to the rear columns 
by rigid elements with all three translational degrees of freedom coupled. However, at 
the from panel, the relative displacement of front panel and front columns are greatly 
dependent on the type of connecting elements used. It should also be noted that the single 
point constraint forces developed at the front columns are also dependent on the amount 
of load that has been transfered to the columns through the front connecting elements. 

Properties of each type of connecting element are chosen in such a way that its de- 
formation can follow the elongated hole displacement curve as close as possible. The gap 
length and frictional force of the contacting surface were chosen arbitrarily at 0.5 inches 
and 15.0 lbs respectively. 

Displacement of the left front interfaces (grid points 17 and 41 were identical and 
were tabulated in Table 1). The single point constraint force of the left front column (grid 
points 1005 and 1006) are tabulated in Table 2. Table 3 shows the displacement of the 
right front interfaces (grid points 20, and 44) while Table 4 shows a single point constraint 
force of the right front column (grid points 1007 and 1008). 
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TABLE 1 DISPLACEMENT VECTOR OF GRID POINT 

LEFT FRONT COLUME 
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TABLE 3 DISPLACEMENT VECTOR OF GRID POINT 
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DISCUSION 


Examination of the displacement of the left and right front panels of the box and the 
load distribution over four supporting columns indicates that types R123 and R23 are the 
two extreme cases. For rigid R123 types, the load was distributed evenly over four columns; 
meanwhile, for rigid R23 types, the load was carried by the rear columns only, and the 
box acts as a cantilever beam. The exact load distribution pattern is somewhere between 
the boundary established by these two types of connections. With finer load increments, 
the displacement curve of a stress dependent material element can be made to match the 
displacement curve of an elongated hole and, therefore, used as a reference to evaluate 
performance of other types of connections. 

In general, static analysis with GAP elements and piecewise linear analysis employ 
an iterative scheme for solution. .They are costly to run especially when models become 
complex. For the problem at hand, the GAP element did not yield a corrected solution 
when under tension load, meanwhile, stress dependent material elements did not perform 
well when under compression. Therefore, these two types of elements were not suitable for 
a model with multiple loading cases where the directions of the load at the interfaces are 
often unknown. 

A couple type connection shows the same characteristics as the R23 type, with the 
displacement shifted due to externally applied couple. It is also difficult to use when the 
direction of the loads at the connection is unknown. However, the solution can be obtained 
with less computing time. 

The solution for rigid elements and spring elements always contains some degree of 
error, but can be obtained with less computing time, and can be used in problems involv- 
ing multiple loading conditions. The spring rate can be adjusted to control the relative 
displacement of the two connecting grid points. 

RECOMMENDATION 

It is recommended that for problems with multiple loading conditions, models should 
be run first using R123 type connections. Then, for interfaces that develop forces larger 
than the frictional forces of the elongated hole, the type R123 connection should be re- 
placed with appropriate spring elements. This methodology yields reasonable results with 
minimum computer run time. 
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Hierarchical Tapered Bar Elements 
Undergoing Axial Deformation 

N. G ones an and S. K. Thampi 
GE Government Services , Houston, Texas 

ABSTRACT: A method is described to model the dynamics of tapered 
axial bars of various cross sections based on the well-known Cralg/Bampton 
component mode synthesis technique. This element is formed in terms of 
the static constraint modes and interface restrained normal modes. This is 
in contrast with the finite elements as implemented in NASTRAN where the 
interface restrained normal modes are neglected. These normal modes are 
in terms of Bessel functions. Restoration of a few of these modes leads to 
higher accuracy with fewer generalized coordinates. The proposed models 
are hierarchical so that all lower order element matrices are embedded in 
higher order element matrices. The advantages of this formulation compared 
to standard NASTRAN truss element formulation Me demonstrated through 
simple numerical examples. 


1. Introduction: Tapered bars and beams have high strength to weight ratios as well 
as architectural advantages. They are frequently employed to model structures in di- 
verse applications, such as ship masts, turbine blades, chimney structures or complex 
frame constructions. NASA (Langley) has tested a truss structure which is made from 
tapered members to be used in space applications [1]. The technical literature on tapered 
beams is indeed vast with a long history [2-10]. Tapered beam finite elements are ei- 
ther simple elements (e.g., Lindberg[2], Rouch/Kao[3]) having two degrees of freedom at 
each end or higher order elements (e.g., Thomas/Dokumaci[4], To[5]) having more than 
four degrees of freedom. Ovunck[6], Avakian/Beskos[7], Gupta[8], Banerjee j Williams [9] 
and Spyrakos/ Chen[10] have used frequency dependent finite elments in their analysis of 
tapered bars. B anerjee/ Williams [9] have developed exact dynamic stiffness matrices for 
Bernoulli- Euler beams. However the approach in References [6-10] involves the unknown 
frequencies of the overall structure. The general framework developed by Engels[ll] and 
applied in References [12-14] allows for the derivation of hierarchical finite elements for 
any type of structural element. This approach does not require the prior knowledge of 
system frequencies, thus overcoming the need for an iterative procedure to compute the 
structural response. In the present paper, a dynamic finite element model for a certain 
class of tapered bars with loads acting only in the axial direction is developed. The ele- 
ment matrices are presented in parametric form and can be easily extended to formulate 
the finite elments for a wide variety of tapered bars covering most practical cases. The 
convergence properties of this dynamic finite element, when compared to the regular finite 
element, are examined using numerical examples. 

2. Ass umed Modes Method: The Lagrangian elastic displacement vector e(x,y,z,t) 
for a generic element can be written as the sum of two separate displacement vectors ej 


124 



and e. 


e = e/ + e 0 ) 

where ej is a quasi-static displacement vector due to the interface displacements q T and is 
expressed as a linear combination of static constraint modes ^/j 

ej = <f>i q! ( 2 ) 


The second part e represents the remainder of the total displacement vector e. It is that 
part of e which is measured relative to e/ by an absolute observer. Clearly, e vanishes 
at the f/ coordinates and therefore can be expressed as a linear combination of assumed 
modes ^ which are restrained at those qj coordinates 

e = 4>q ( 3 ) 


The vector q represents a_set of generalized coordinates to be determined as part of the 
solution. One example of ^ modes are the normal modes of the element E restrained at the 
qj coordinates. It should be stressed that although restrained normal modes have often 
unique advantages, they are only one of many possible sets. In fact, (j> modes need only be 
restricted to admissible functions that vanish at the q t coordinates. Substituting Eqs. (2) 
and (3) into Eq. (1) yields 


e = 




( 4 ) 


l 9 


so that e is written in terms of a linear combination of two sets of assumed modes: (1) 
static constraint modes and (2) interface restrained assumed modes. It should be noted 
that the representation of e in Eq. (4) is complete in the sense that any degree of accuracy 
is theoretically possible as long as enough <f> modes are added. 

In the standard consistent mass matrix approach, the elastic displacement vector e 
over the element is represented as a linear combination of interpolation or shape functions. 
In fact, these shape functions are identical to the static constraint modes <j>i and therefore e 
is approximated by «/ as in Eq. (1). The standard finite element approach therefore totally 
neglects the displacement e in Eq. (1). Ignoring e leads directly to a deterioration of the 
modal content of a typical finite element model. One way to ensure better convergence 
to a desired model fidelity is suggested by Eq. (4) and leads to dynamic finite element 
models. Instead of totally neglecting the e displacement, one could retain a limited number 
of q coordinates, thereby improving the mass and force distribution models. Of course 
adding q coordinates also increases the order of the overall model. However, this approach 
has three important advantages: (1) The model converges much faster, i.e., far fewer 
degrees of freedom are necessary to attain comparable accuracy; (2) in principle, no further 
subdivision of basic elements is necessary, thereby simplifying the finite element grid and 
(3) the model is hierarchical and therefore has all the advantages associated with this 
property. In addition, these finite element models are directly based on the assumed 
modes method which provides a sound theoretical basis. 
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In Reference [11] it is shown that for a linear elastic 
potential energies T and V can be written as 


material, the element kinetic and 


T = 




(5) 


V = \l T Xq, K ~ j (B4,) T C(B4>) dV (6) 

E 

in which the matrix C is the material stiffness matrix and B contains the appropriate 
partial derivatives in x,y and z. Futhermore, the matrices M and K represent the mass 
and stiffness matrices of a generic finite element E. In partitioned form, 


Mji 

Min 

K = 

Ku 

0 ' 

-MJn 

Af/vivJ ’ 

0 

Kni v. 


(7) 


where 


and 


Mu = J <t>J 4>i dm, Min = J 4>J<t> dm, Mnn = J 

E E E 


J 

Af/y/v = / <f> 4> dm 


Kii — J (B<j>i) T C(B<f>i)dV, K NN = J(B<t>) T C(B<f>)dV 


( 8 ) 

(9) 


Note that the Kin partition is always zero, which means that no stiffness coupling exists 
between qi and q. 


At this point, a few remarks are in order. First, the matrices Mu and Ku repre- 
sent the standard finite element consistent mass and stiffness matrices for the element E. 
The consistent mass matrix approach represents in fact a static condensation or Guyan 
reduction whereby all noninterface degrees of freedom are eliminated. Secondly, if the in- 
terface restrained normal modes are used for the columns of <j>, then the present approach 
is identical to the Craig/Bampton component mode synthesis procedure as applied to a 
finite element. It should be emphasized that the element E is generic. This means that 
the proposed approach is valid, at least in theory, for any type of element. In the present 
paper, this general procedure will be applied to the special case of the tapered bars. 


3. Tapered Bars: The hierarchical stiffness and mass matrices of the tapered bar are 
obtained by solving the governing equations of motion for displacement. Figure (1) repre- 
sents linear tapered bar ab with a straight centroidal axis and the directions of the principal 
axes being the same for all crossections. The cross sectional area A(x) is given by 

A{x) = A a (l +cjY (10) 

where c = db/d a — 1 and A,-, <f, (» = a,b ) denote the cross sectional area and the depth 
respectively, c > — 1 otherwise the beam tapers to zero between its ends and L is the 
length of the bar. 
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Although the formulation is valid for any n > 0, many 
practical cases of tapered bars arise when n is one or two, see 
Figure (2). If the geometrical properties of the element at both 
ends are given, the shape function for n can be derived as 

_ log{AbJ_^a) ( 11 ) 

n ~ log{d b /d a ) 

For bars of closed box or I-section of constant width and vary- 
ing depth, n is not an integer and will vary slightly from E- 
q. (10) at all * other than the two ends. But the deviations 
are usually within one percent of the exact values. 



Tapered Element 


4. Static Constraint Modes: Consider the axial bar element as illustrated in Fig- 
ure (1). The bar is assumed to undergo vibration along its own axis and as a rigid o y can 
only move along that same axis. The interface displacements are defined as 9l (t) - 
and q 2 (t) = u(L,t), i.e., q t = [ «*i ? 2 f and the displacement vector e is considered to have 
only one component, i.e., u. Eq. (4) is therefore written as 


u — ( <h <f >2 Jfl/ + 4 9 


( 12 ) 



(a) <D) (c) <d) (e) 



Figure 2. Sample Cross Sections 
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The static constraint modes, 4>\ can be computed from the differential equation for 
the axial deformation u at * from end a of the tapered bar, 




where E is Young’s modulus. Integration of Eq. (13) gives, 

_ .. .du „ 

EA Wj; = Ci 

Substituting for A(x) from Eq. (10) and setting ( = 1 + c jr gives 

du Ci L 1 


Integrating Eq. (15) again, 


where 


u — 


di EA a c ( n 

CiL 


EA a 


-m+c 2 


(13) 


(14) 


(15) 


(16) 


m = -■ 


— for n ^ 1 


(17) 


(n - 1)** 

= In £ for n = 1 

The appropriate boundary conditions for the computation of static constraint modes are 
given by 

,l = l ’ fc= ° (18) 


9i = 0, <72 = 1 

The resulting static constraint modes are 


^ /(0-/(i + c) /(D-/(0 


(19) 


Note that these constraint modes are in fact the shape functions used in the stiffness matrix 
of the tapered bar. 

At this point, enough information exists to compute Mjj and Ku from Eqs. (8-9). 
Indeed, for the tapered bar, the kinetic and potential energies T and V are given as 




( 20 ) 


Substituting Eq. (12) into Eq. (20) gives 

L 


Mu = j pA(x)<f>J(f>idx, Ku = Jem x) 


d<f>i T d<f>i 
dx dx 


dx 


( 21 ) 
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When n = 1, the explicit expressions for Ku and Mu are defined as 


Kn = 


EA a c 
L ln(l + c) 



-1 

1 


Mu = 


m ii 

m-n 


m 12 
m 22 


(22) 


m„ = ft 1 ~ s t- 21 "^ 1 + <0 - 21 ”(1 + <0 + + 2c l 

4cln (1 + c) 

mi 2 = m 2 i = - T [(c 2 +2C+ 2)In(l + c) - c 2 - 2,] (23) 

4cln (1 + c) 

m ,2 = . [2(1 + c) 2 (ln 2 (l + c) - ln(l + c)) + c 2 + 2c] 

4cln (1 + c) 

The counterpart expression for Ku when n ^ 1 is given 


Ku = 


EA a c(n - 1)(1 + c) 


n — 1 


1 -l 

-l l 


(n^l) 


(1 + C )"- 1 - 1 

The stifness matrix partitions Kji and Mu for n = 2 are evaluated as 

„ ^A a (l + c) 

Ku = 


(24) 


1 -1* 

I __ pA a L 1 

2 

1 + C 

-1 1 

& 

II 

as 

1 +c 

2(1 + c) 2 


(25) 


Because of the similarity of the governing equations between axial bars and torsional 
shafts, Eq. (25) can be used as the stiffness and consistent mass matrices for tapered shafts 
by replacing the variables A a , E with J a , G respectively, where J a is the polar second 
moment of area and G is the shear modulus of elasticity. If it is decided that no extra q 
coordinates are to be retained in Eq. (4), then the procedure can be terminated at this 

stage. 

S. Interface Restrained Normal Modes: An entire class of hierarchical models can 
now be created solely on the basis of choosing the interface restrained assumed modes. In 
this paper, the set of interface restrained normal modes is used. The normal modes and 
their corresponding frequencies are obtained from solving the eigenvalue problem associ- 
ated with the partial differential equation, 


d_ 

dx 



(26) 


subjected to the clamped-clamped boundary conditions. 

Substituting for A(x) from Eq. (10) and letting ( = l + cx/L gives, 

Pu ndu pL 2 0 2 u 
8 + Tdt " Ec 2 di 2 


(27) 
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For harmonic vibration, 


u(£,*) = {/(£)sinu>< 


(28) 


where t denotes time and u> is the circular natural frequency. Eq. (27) is modified using 
Eq. (28) as 

cPU ndU 


+ = 0 


Ec 2 


(29) 


The solution of Eq. (29) when c > 0 is 

U(Q = (^ + 


(30) 


where J and Y are Bessel functions of the first and second kind and a = — « Pk. 

c v E 

For the case of n = 2, imposing the clamped-clamped boundary conditions in Eq. (30) 
yields the characterstic equation for the tapered bar, 


sin ac = 0 


The solution of Eq. (31) is 


(31) 


» 1 = 1,2 oo (32) 

It is noteworthy that the natural frequencies of the tapered bar with clamped ends for the 
case of n = 2 are independent of c and are in fact the same as that of uniform rods. From 
Eq. (30), the interface restrained normal modes for n = 2 becomes 



= C ‘ i"K~°) 

where the mass normalization constant (7, must satisfy 

1 + C 

/ = i 

1 

Substituting Eq. (33) into Eq. (34), C{ is determined as 

c ‘=v Gz 

From Eq. (8), the mass matrix partitions are evaluated as 

L l 

Mnn = J pA{x]4> J <f>dx - I , Mi N = [mj,] , m,ji - J pA(x)<f> j ^ i dx 

o o 

i 1,2,..., oo , j = 1)2 


(33) 


(34) 


(35) 


(36) 
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with 


mi,- = 


■y/2 pA a L 


m 2 i = 




( 1 + c)(-iy 


When c is set to zero, Eqs. (25) and (36) reduce to the case of uniform bars. Appli- 
cation of Eq. (33) into Eq. (20) gives the stiffness matrix partitions, 

Kin = 0, Knn — Dmg* [w* . . .] (38) 


and ui is given by Eq. (32). It can be seen that the tapered bar element matrices consist 
of very simple terms. 


6. Demonstration Examples: 


6.1 Tapered Cantilever: This first example is concerned 
with the cantilevered bar clamped at the right end, see Fig- 
ure (3). The pertinent structural parameters are E = 30 x 
10® psi, pg = 0.2839 lb/in 3 ,A a = 1 in 2 ,A b = 4 in 2 and L = 
72 in. 

The characteristic equation of this cantilever is given by 



a cos ac + sin ac = 0 (39) 

Figure 3. Cantilever 

Tables (1) and (2) show the number of converged frequencies 
to within a given percentage relative error when compared to 
the theoretical frequencies from Eq. (39) for different model 
orders n. 

Note that the bar is subdivided into the requisite number of elements to arrive at the 
model order n in the standard finite element method as implemented in CSA-NASTRAN 
whereas in the hierarchical finite element model, the number of bar elements is always one 
and the requisite interface restrained normal modes are added to arrive at n. 


Table 1. Standard Finite Element Method 
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Table 2. Hierarchical Finite Element Method 


e 

71 

< i % 

< 5 % 

< 10 % 

8 

16 

24 

32 

40 


14 

23 

31 

39 

m 

15 



39 

D 

15 

23 

31 

39 


6.2 Planar Truss: Next, consider the planar single-bay truss 
as illustrated in Figure (4). The parameters are the same as 
in Section (6.1). The horizontal bar has a uniform cross sec- 
tional area of 4 in 2 . For planar truss elements, the transverse 
inertia must be taken into account. This time, however, no 
‘exact solution’ for the frequencies of the structure exists. A 
reference solution was obtained in two different ways: (1) by 
constructing a highly refined standard finite element model (2) 
by retaining a large number of normal modes in the hierarchical 
model. Both models were refined to the point where no signif- 
icant change in the frequencies occured and both the models 
produced the same results. 



Figure 4. Planar Truss 


Table (3) lists the frequencies from the hierarchical finite element model when two 
normal modes per bar are added (i.e., n = 16) as compared to the reference solution. In 
order to achieve comparable results from standard finite element method, each bar has to 
be subdivided into five beam elements (i.e., n = 66) and consistent mass matrix has to be 
generated within NASTRAN. 


Table 3. Frequency Comparisons For Planar Truss 


Mode 

Number 

Reference 
Frequency (Hz) 

Computed 
Frequency (Hz) 

Error 

% 

2 

2.606E2 

2.606E2 

-1.123E-2 

3 

2.999E2 

2.999E2 

-1.703E-2 

5 

6.958E2 

6.961E2 

-4.763E-2 

7 

1.467E3 

1.468E3 

-1.837E-2 

9 

1.757E3 

1.762E3 

-2.970E-1 

11 

2.141E3 

2.172E3 

-1.432 

13 

2.889E3 

2.898E3 

-3.317E-1 

14 

3.049E3 

3.081E3 

-1.045 

1 16 

3.382E3 

3.505E3 

-3.622 
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7. Conclusions For the case of tapered bars, it has been demonstated that the modal 
synthesis approach, where substructures are assembled to form the overall structure, can 
be used at the element level itself. This approach has superior convergence characteristics 
and the advantages of hierarchical formulation. 
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ABSTRACT 

A method for computing transient thermal stress vectors from 
temperature vectors is described. The three step procedure 
involves the use of NASTRAN to generate an influence coefficient 
matrix which relates temperatures to stresses in the structural 
model. The transient thermal stresses are then recovered and 
sorted for maximum and minimum values. Verification data for 
the procedure is also provided. 


1-0 INTRODUCTION 

There are many occasions when the stresses produced by transient 
thermal events must be considered. The ascent, on-orbit, and 
descent phases of a Spacelab mission produce large temperature 
< 3*'<^dients on the Cargo Element. A method for recovering the 
thermal stresses produced by these events was developed by the 
Structural Analysis Group at McDonnell Douglas Space Systems 
Company-Hunts v i 1 le (MDSSC-HSV) , and has been used for more than 
six years in Spacelab Evaluation. 

Because this method was somewhat unhandy to use, only a limited 
number of temperature distributions could be run. These were 
generally chosen on the basis of temperature or temperature 
difference extremes. Unfortunately, there is no proven 
relationship between "worst stresses" and "worst temperatures" 
or "worst temperature differences" for the complex models with 
which Spacelab has to deal. It was therefore decided to develop 
a process that would transform transient temperature vectors 
(which are separately calculated) directly into stress vectors 
that could then be sorted for maximum and minimum (Max/Min) 
values. 

Section 2 presents the theory used in the procedure, Section 3 
describes the procedure, Section 4 presents verification 
results, and Section 5 contains the conclusions. 
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2.0 


THEORY 


The method described assumes that the thermal strains, and hence 
stresses are linearly proportional to the temperatures. The 
method, therefore, is not applicable to models having 
temperature or stress dependent material properties. 

An influence coefficient matrix, [TRAS], to transform the 
temperatures at a models n grids into m stresses can be 
developed as follows. A temperature of one unit above the 
reference is applied to a grid in the model while the remaining 
n-1 grids are held at the reference temperature. The resulting 
m stresses are then recovered. This is repeated for all n grids 
in the model. The resulting n sets of stresses are then 
assembled as columns to form an m by n influence coefficient 
matrix, [TRAS], that can be used to transform temperature 
vectors into stress vectors. The transformation is as follows: 


Where ; 


[ STMHST ] = [TRAS] [TEMPS] 


( 1 ) 


[STMHST] = Stresses in the finite element model 
(time history) . 

[TRAS] = Linear transformation (influence 
coefficient) from temperatures to 
stresses . 

[TEMPS] = Temperatures at the grids (time 
history) . 


Equation (1) is used to recover transient stress vectors from 
transient temperature vectors. 


3 . o PROCEDURE 

The procedure used to recover maximum and minimum thermal 
stresses is divided into three steps. The first step is to 
generate the influence coefficient matrix [TRAS] . The transient 
stresses are recovered using equation (1) in the second step. 

The transient thermal stresses are then sorted for Max/Min data 
in the third step. Each of these three steps will be described 
below. 
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3.1 


Step 1 - Generation of The Influence 
Coefficient Matrix [TRAS] 

NASTRAN is used to generate [TRAS] as described in Section 2. 
NASTRAN Subcase commands are used to accomplish this step with 
each subcase corresponding to a grid in the model. It should be 
noted that NASTRAN can handle a maximum of sixty-six temperature 
load cases. A model having more than sixty-six grids will, 
therefore, require multiple runs and subsequent merging of the 
output data. An example of a NASTRAN control deck to generate 
[TRAS] for the verification model described in Section 4 is 
shown in Figure 1. The DUMM0D5 module is used to convert the 
OES1 table into a matrix data block which is then written to a 
file using the OUTPUTS module. The extraneous rows (fiber 
distances, safety margins, etc.) of [TRAS] are then removed 
using a specially developed FORTRAN code. 


3-2 Step 2 - Recover Transient Thermal Stress 

Vectors From Temperature Vectors 

The transient temperature vectors must now be obtained. At 
MDSSC-Huntsville, the thermal analyses are performed using 
SINDA. A FORTRAN code has been written to access the SINDA 
output file and recover the desired temperatures along with the 
corresponding time vectors. The temperature and time data is 
written in OUTPUTS format. Care should be taken to ensure that 
the row order of the temperature vectors is compatible with the 
column order of [TRAS]. 

Equation (1) is then used to recover the transient thermal 
stress vectors. This is easily accomplished in a simple DMAP 
run using the MPYAD module. The resulting thermal stress 
vectors and time vector are written to a file using the OUTPUTS 
module. 


3-3 Step 3 - Sort Thermal Stress Vectors For 

Maximum and Minimum Data. 

The transient thermal stress vectors are now sorted for Max/Min 
data. This step is performed using a specialized FORTRAN code. 
This code can search multiple time histories allowing composite 
Max/Min tables to be obtained. An example of output from this 
program for the Spacelab Multipurpose Experiment Support 
Structure (MPESS) is shown in Table 1. The data in Table 1 was 
obtained from actual temperature vectors for a Spacelab mission 
and encompasses ascent, orbit, and descent. 
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4.0 


VERIFICATION 


In order to verify that the procedure is working correctly, a 
test case was performed. A simple rectangular plate model was 
constructed using QUAD4 and BAR elements. A plot of the plate 
model is presented in Figure 2 and the Bulk Data is shown in 
Figure 3. The model has thirty-three grids and is homogeneous. 

A temperature gradient (see Figure 4) was applied to the model 
and the resulting stresses in four elements were recovered using 
NASTRAN directly. The results are presented in Figure 5. The 
transient thermal stress recovery procedure was then used to 
calculate stresses due to the same temperature gradient and the 
results are shown in Table 2 . 

It can be seen from Figure 5 and Table 2 that the results are 
the same. This indicates that using [TRAS] to perform the 
thermal stress recovery produces the same results as using 
NASTRAN directly. 


5.0 CONCLUSION 

A procedure has been developed to recover thermal stresses in a 
NASTRAN model directly from temperatures output by a thermal 
model. Because the procedure uses a linear transformation 
matrix rather than a computer program, entire thermal stress 
time histories may be efficiently obtained and scanned for 
Max/Min thermal stress data. Tabular output of the Max/Min data 
is then produced. A test case has been executed and the results 
indicate that the procedure functions correctly. 
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NASTRAN TITLEOPT=0 

$ 

$ RUN TO GENERATE INFLUENCE COEFF. MATRIX FOR THERMAL STRESS RECOVERY. 

$ 

$ 

$ 

S 

ID CHECKOUT PLATE 
APP DISP 
SOL 1,7 
TIME 999 

DIAG 13, 14,21 ,22,26,42 
$ 

$ WRITE OUT ELEMENT STRESS/FORCE MATRICES 
$ WHERE; 

$ OESl=STRESSES 

S OEFl=FORCES 

$ I =# OF ELEMENTS FOR WHICH STRESSES/FORCES ARE BEING RECOVERED 

$ 

ALTER 106 

DUMM0D5 OES1 , , , ,/ELSEU, , , ,/C,N,AA/////C,N, 1/C,N, 1 $ 

OUTPUT5 ELSEU, , , , //0/15//0 $ 

OUTPUT5 , , , .//-9/15//0 $ 

ENDALTER 

CEND 

$ 

TITLE = ANALYSIS OF PLATE 

SUBTITLE - RUN TO GENERATE INFLUENCE COEFFICIENT MATRIX 
LABEL = PLATE 
MAXLINES = 99999999 
SPC=100 

$ 

SUBCASE 1 

LABEL = APPLY A 1 DEGREE ABOVE REF. TEMP. TO GRID 1 

TEMP (LOAD) = 1 

ELSTRESS=ALL 

$ 

SUBCASE 2 

LABEL = APPLY A 1 DEGREE ABOVE REF. TEMP. TO GRID 2 

TEMP (LOAD) = 2 

ELSTRESS=ALL 

$ 

SUBCASE 3 

LABEL = APPLY A 1 DEGREE ABOVE REF. TEMP. TO GRID 3 

TEMP (LOAD) = 3 

ELSTRESS=ALL 

$ 

SUBCASE A 

LABEL = APPLY A 1 DEGREE ABOVE REF. TEMP. TO GRID A 

TEMP (LOAD) = A 

ELSTRESS=ALL 

$ 

SUBCASE 5 

LABEL - APPLY A 1 DEGREE ABOVE REF. TEMP. TO GRID 5 
Figure 1. NASTRAN Control Deck to Generate [TRAS] for the 
Verification Plate Model 
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TEMP (LOAD) - 5 
ELSTRESS-ALL 

$ 

SUBCASE 6 

LABEL - APPLY A 1 DEGREE ABOVE REF. TEMP. TO GRID 6 

TEMP (LOAD) - 6 

ELSTRESS-ALL 

$ 

SUBCASE 7 

LABEL = APPLY A 1 DEGREE ABOVE REF. TEMP. TO GRID 7 

TEMP (LOAD) - 7 

ELSTRESS-ALL 

$ 

SUBCASE 8 

LABEL = APPLY A 1 DEGREE ABOVE REF. TEMP. TO GRID 8 

TEMP (LOAD) - 8 

ELSTRESS-ALL 

$ 

SUBCASE 9 

LABEL - APPLY A 1 DEGREE ABOVE REF. TEMP. TO GRID 9 

TEMP (LOAD) - 9 

ELSTRESS-ALL 


$ 

SUBCASE 10 

LABEL = APPLY A 1 DEGREE ABOVE REF. TEMP. TO GRID 10 

TEMP (LOAD) = 10 

ELSTRESS-ALL 

$ 

SUBCASE 1 1 

LABEL - APPLY A 1 DEGREE ABOVE REF. TEMP. TO GRID 11 

TEMP (LOAD) - 11 

ELSTRESS-ALL 

$ 

SUBCASE 12 

LABEL - APPLY A 1 DEGREE ABOVE REF. TEMP. TO GRID 12 

TEMP (LOAD) - 12 

ELSTRESS-ALL 

$ 

SUBCASE 13 

LABEL - APPLY A 1 DEGREE ABOVE REF. TEMP. TO GRID 13 

TEMP (LOAD) - 13 

ELSTRESS-ALL 

$ 

SUBCASE 14 

LABEL - APPLY A 1 DEGREE ABOVE REF. TEMP. TO GRID 14 

TEMP (LOAD) - 14 

ELSTRESS-ALL 

$ 

SUBCASE 15 

LABEL - APPLY A 1 DEGREE ABOVE REF. TEMP. TO GRID 15 

TEMP (LOAD) - 15 

ELSTRESS-ALL 


SUBCASE 16 

LABEL = APPLY A 1 DEGREE ABOVE REF. TpiP. T0 . £ RI Dl|c i for the 
Figure 1. NASTRAN Control Deck to Generate [TRAS] 
Verification Plate Model (Continued) 
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TEMP (LOAD) = 16 
ELSTRESS=ALL 


$ 

SUBCASE 17 

LABEL - APPLY A 1 DEGREE ABOVE REF. TEMP. TO GRID 17 

TEMP (LOAD) = 17 

ELSTRESS=ALL 

$ 

SUBCASE 18 

LABEL - APPLY A 1 DEGREE ABOVE REF. TEMP. TO GRID 18 

TEMP (LOAD) = 18 

ELSTRESS=ALL 

$ 

SUBCASE 19 

LABEL - APPLY A 1 DEGREE ABOVE REF. TEMP. TO GRID 19 

TEMP (LOAD) = 19 

ELSTRESS-ALL 

$ 

SUBCASE 20 

LABEL = APPLY A 1 DEGREE ABOVE REF. TEMP. TO GRID 20 

TEMP (LOAD) = 20 

ELSTRESS=ALL 

$ 

SUBCASE 21 

LABEL = APPLY A 1 DEGREE ABOVE REF. TEMP. TO GRID 21 

TEMP (LOAD) = 21 

ELSTRESS=ALL 

$ 

SUBCASE 22 

LABEL - APPLY A 1 DEGREE ABOVE REF. TEMP. TO GRID 22 

TEMP (LOAD) =22 

ELSTRESS=ALL 

$ 

SUBCASE 23 

LABEL - APPLY A 1 DEGREE ABOVE REF. TEMP. TO GRID 23 

TEMP (LOAD) = 23 

ELSTRESS-ALL 

$ 

SUBCASE 24 

LABEL = APPLY A 1 DEGREE ABOVE REF. TEMP. TO GRID 24 
TEMP (LOAD) = 24 
ELSTRESS=ALL 
5 

SUBCASE 25 

LABEL = APPLY A 1 DEGREE ABOVE REF. TEMP. TO GRID 25 

TEMP (LOAD) = 25 

ELSTRESS-ALL 

$ 

SUBCASE 26 

LABEL - APPLY A 1 DEGREE ABOVE REF. TEMP. TO GRID 26 

TEMP (LOAD) = 26 

ELSTRESS-ALL 

$ 

SUBCASE 27 

LABEL - APPLY A 1 DEGREE ABOVE REF. TEMP. TO GRID 27 

Figure 1. NASTRAN Control Deck to Generate [TRAS] for the 
Verification Plate Model (Continued) 
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TEMP (LOAD) “ 27 
ELSTRESS=ALL 


SUBCASE 28 

LABEL = APPLY A 
TEMP (LOAD) = 28 
ELSTRESS=ALL 

$ 

SUBCASE 29 

LABEL = APPLY A 
TEMP (LOAD) = 29 
ELSTRESS=ALL 

$ 

SUBCASE 30 

LABEL = APPLY A 
TEMP (LOAD) = 30 
ELSTRESS=ALL 

$ 

SUBCASE 31 

LABEL = APPLY A 
TEMP (LOAD) = 31 
ELSTRESS=ALL 

$ 

SUBCASE 32 

LABEL = APPLY A 
TEMP (LOAD) = 32 
ELSTRESS=ALL 

$ 

SUBCASE 33 


LABEL 

= APPLY 

A 1 DEGREE AB 

TEMP (LOAD) « 
ELSTRESS-ALL 

$ 

BEGIN BULK 
$ 

33 


SPC1 

$ 

100 

123456 

1 

TEMP 

1 

1 

71.0 

TEMPD 

1 

70.0 


TEMP 

2 

2 

71.0 

TEMPD 

2 

70.0 


TEMP 

3 

3 

71.0 

TEMPD 

3 

70.0 


TEMP 

4 

4 

71.0 

TEMPD 

4 

70.0 


TEMP 

5 

5 

71.0 

TEMPD 

5 
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TEMP 

6 

6 

71.0 

TEMPD 

6 
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TEMP 

7 

7 

71.0 

TEMPD 

7 

70.0 


TEMP 

8 

8 

71.0 

TEMPD 

8 

70.0 


TEMP 

9 

9 

71.0 

TEMPD 

9 

70.0 


Figure 1 . 

NASTRAN Control 

Verification Pi; 


1 DEGREE ABOVE REF. TEMP. TO GRID 28 


1 DEGREE ABOVE REF. TEMP. TO GRID 29 


1 DEGREE ABOVE REF. TEMP. TO GRID 30 


1 DEGREE ABOVE REF. TEMP. TO GRID 31 


1 DEGREE ABOVE REF. TEMP. TO GRID 32 


REF. TEMP. TO GRID 33 
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Figure 1. NASTRAN Control Deck to Generate [TRAS] for the 
Verification Plate Model (Concluded) 
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Figure 3. Plate Model Bulk Data (Concluded) 
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Figure 5. Thermal Stresses from a Direct NASTRAN Solution of 
the Plate Model 
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Figure 7. MPESS BAR 106 (SINDA Node 6000) Axial Stress Time 
History 
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Figure 8. MPESS BAR 106 (SINDA Node 6000) Temperature Time 
History 
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m>SSC - HUNTSVILLE OPERATIONS 


TABLE 2. THERMAL STRESSES IN THE PLATE MOOEL FROM THE TRANSIENT THERMAL STRESS RECOVERY PROCEDURE 
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Transient Loads Analysis 
For Space Flight Applications 

S, K . Thampij N. 5. Vidyasagar and N. Ganesan 
GE Government Services , Houston , Texas 

ABSTRACT: A significant part of the flight readiness verification process 
involves transient analysis of the coupled Shut tie- payload system to determine 
the low frequency transient loads. This paper describes a methodology for 
transient loads analysis and its implementation for the Spacelab Life Sciences 
Mission. The analysis is carried out using two major software tools - NAS- 
TRAN and an external FORTRAN code called EZTRAN *. This approach is 
adopted to overcome some of the limitations of NASTRAN’s standard tran- 
sient analysis capabilities. The method uses Data Recovery Matrices (DRM) 
to improve computational efficiency. The mode acceleration method is ful- 
ly implemented in the DRM formulation to recover accurate displacements, 
stresses and forces. The advantages of the method are demonstrated through 
a numerical example, 

~ — Introduction: In the past decade, NASA has conducted numerous Spacelab Missions 
for the advancement of space exploration and research. The Spacelab is a reusable labora- 
tory that is carried in the cargo bay of the Space Shuttle Orbiter. Experiments in several 
different disciplines such as astronomy, life sciences and material science are accommodated 
in this modular laboratory for various Shuttle missions. The module also contains utili- 
ties, computers and work benches to support the experiments. The experiment hardware 
is mounted in instrument racks located on either side of the module, in overhead lockers, 
and in the center aisle, as shown in Figure 1. 

During liftoff and landing flight events, the Shuttle and its payload show significant 
low-frequency transient accelerations due to thrust from the main engines and solid rocket 
boosters, wind gust, vortex shedding, and launch pad forces during liftoff, and crosswinds 
and nose-gear slapdown during landing. The levels of acceleration on a specific payload 
component depend on the response of the Spacelab inside the Orbiter cargo bay and 
the response of the component inside the Spacelab. Because these responses depend on 
the dynamic characteristics and interactions of the Orbiter-Spacelab-payload system, a 
transient analysis of the coupled system is required to determine the quasi-static loads as 
part of the flight readiness verification process. 

Analysis of the coupled system can be carried out using the standard transient analysis 
capabilities of NASTRAN. However, these procedures have some limitations in terms of 
computational efficiency and accuracy, especially when dealing with large substructured 
models [Rf. 1,2]. An alternate procedure which relies on the extensive use of data recovery 
matrices is presented to overcome some of these limitations. The methodology has been 
successfully implemented for the analysis of the Spacelab Life Sciences (SLS-2) Mission [Rf. 


* EZTRAN is developed by Structural Dynamics Research Corporation, San Diego, 
California 
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3] which involved large models (in excess of 190000 degrees of freedom) and long simulations 
(in excess of 4500 time steps). The advantages of the method are demonstrated through a 
simple example problem in this paper. 

2. Overview of Methods: There are two methods for performing dynamic transient 
analysis in NASTRAN. The Direct Transient Response, available in Rigid Format 9, solves 
a system described in terms of its physical mass, stiffness and damping matrices. 

[m] {u} + [c] {«} + [ k ] {«} = {/} (1) 

These equations are numerically integrated to determine the response at the physical de- 
grees of freedom (DOF) as functions of time. The Modal Transient Response method, 
available in Rigid Format 12, is similar to the direct method, with the exception that it 
uses the classical modal transformation 

{«} = [<f>] {Q} ( 2 ) 

to diagonalize the physical mass, stiffness and damping matrices. The overall response is 
calculated by including only a small number of the structural modes, making the numerical 
integration of the generalized equations of motion much faster. 

IM] {«} +[<?]{«} + [*]{«} = { f > < 3 > 

where 

[M] = [<i>] T [m] [<f >] , [C] = [4>f [c] M , [K] = [<f>] T [*] [<j>] , {F} = [<t>\ T {/> 

Physical responses can be recovered from a modal transient analysis through the mode 
displacement method (Eqn. 2) or the mode acceleration method (Eqn. 4). The former 
method is more efficient and quite accurate for calculating accelerations if sufficient modes 
are retained to envelope the frequency content of the forcing functions. However, the 
latter method is preferred when accurate displacements, element forces and stresses are 
required from a modal transient analysis. The mode acceleration technique minimizes the 
loss of accuracy due to modal truncation by including the static response of the truncated 
high-frequency modes in the solution. 

M = [fe ]- 1 [{/} - M {*} - M {«}] ( 4 ) 

The implementation of the mode acceleration procedure in NASTRAN has some disad- 
vantages. In order to include the effects of inertia and damping forces in the load vector 
in Eqn. 4, the accelerations and velocities at all DOF in the solution set (L-set) must be 
computed. A static solution must then be performed with the modified load vector at each 
time step. This requires significant computer processing time if the L-set is large and/or 
the number of integration steps is large. In addition, there are accuracy problems when 
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dealing with multi-stage substructured models. For such cases, the mode- acceleration cor- 
rection is applied only to the residual structure and not to upstream substructures. This 
leads to modal truncation errors because the component modes are used to represent the 
static response of the interior degrees of freedom [Rf. 1]. 

Because of these disadvantages, standard transient analysis procedures in NASTRAN 
are not suitable for solving large problems involving multi-stage substructured models. 
Alternate methods are required to overcome these limitations, and yet retain the benefits 
of the mode acceleration method. Such a procedure, used for transient loads analysis of 
the SLS-2 Mission configuration, is described in the following section. 

3. Alternate Method: The alternate method is based on a slightly different form of the 
mode acceleration data recovery equation. Assuming that damping is negligible, Eqn. 4 
can be expressed in the following convenient form. 

{«} = M {p} - M [A -1 ] {<$} (5) 

The first term represents the static portion of the transient response. It is obtained as the 
product of the static response caused by a set of unit loads, I, and a time varying load 
scale factor, p(f). Note that the unit loads multiplied by the scale factor are equivalent 
to the applied loads, i.e., 

{/(<)} = M W0> (6) 

The second term in Eqn. 5 represents the dynamic portion of the response. It usually 
includes only the elastic mode contributions as the rigid body modes (if any) do not 
contribute to stresses and forces. However, the contribution from rigid body modes to the 
total displacement, [<£,.(,] {Q r t}, can be included in Eqn. 5 if desired. The computational 
advantages of the alternate method stem from the size of the [V>] and [<j>] matrices which 
are determined by the number of response recovery points, the number of load application 
points and the number of retained elastic modes. These are usually much smaller than the 
full model size. 

The alternate method is implemented using two major software tools - NASTRAN 
and an external FORTRAN program called EZTRAN. NASTRAN is used to develop, 
process and assemble the finite element model of the coupled system, calculate system 
modes, determine unit load static responses, and create data recovery matrices. EZTRAN 
calculates the modal initial conditions, solves the generalized equations of motion, and 
recovers physical results. The following steps describe how the two work in conjunction to 
perform the various analysis tasks. 

3.1 Model Generation and Assembly: Finite element models of the Orbiter, S- 

pacelab and experiment payloads are developed by different organizations and are usually 
test-verified models. They are assembled into a solution system using the automated 
multi-6tage substructuring features of NASTRAN. Prior to assembly, the quality of each 
component model is verified by performing a series of analytical checks including rigid 
body modes check, stiffness matrix equilibrium check, rigid body mass check and an en- 
forced displacement check. At each stage of assembly, the effective DOF in the model is 
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reduced either through a Guyan reduction (REDUCE) or a modal reduction (MREDUCE). 
A fixed-interface Craig-Bampton modal reduction is the preferred method as it will not 
compromise the fidelity of test-verified finite element models. 

3.2 Loads Definition: The input forcing functions for coupled loads analysis are ob- 
tained from previous flight accelerometer data. These are maintained by NASA and pro- 
vided to payload organizations during design evaluation to refine design loads. The data is 
provided in terms of discrete force-time coordinate pairs, with the point of application of 
each force being identified by a node number and component name. They include multiple 
load cases for liftoff and landing flight events. Load scale factors are generated for each 
load case by normalizing the system forcing functions with unit loads, as indicated in Eqn. 
6. The unit loads are defined as the maximum force occuring at each loaded DOF, across 
all load cases. 

3.3 Normal Modes Analysis: A normal modes analysis of the fully assembled system is 
performed using NASTRAN Rigid Format 3. The rigid body modes and the elastic modes 
of the system, in a specified frequency range, are recovered and stored in the substructure 
operating file (SOF) database. The frequency range for modal truncation is decided based 
on the frequency content of the excitation. 

3.4 Unit Load Static Analysis: An inertia relief static analysis is performed on the 
fully assembled system for unit loads derived from the liftoff and landing forcing functions. 
An unit load vector is generated for each loaded DOF, and they are sequenced in the same 
order as the forcing functions to form a unit load matrix. The static analysis is performed 
using Rigid Format 2 because it is capable of analyzing structures with rigid body modes. 
SUPORT cards must be included if rigid body modes are present, and the choice of support 
points has significant effect on the computation of displacement results. A good choice 
is indicated by low strain energy at the support points. An unreduced model is used for 
static analysis in order that the full mass matrix be available for calculating internal inertia 
loads of upstream substructures. The static displacements from the inertia relief solution 
are recovered and stored in the SOF database. 

3.5 Data Recovery Matrices Generation: The alternate procedure requires the 
generation of acceleration and displacement data recovery matrices. These are formed for 
each basic substructure by performing two data recovery (Phase 3) restart runs with special 
DMAP alters. The acceleration DRM is made up of rigid body modes and retained elastic 
mode vectors. These are extracted from the normal modes database for DOF specified 
through XYPLOT/XYPEAK requests in the SOL 3 data recovery run. 

Acceleration DRM = [4> r b 4>ei] CO 

Displacements and displacement dependent responses such as element forces, stresses and 
substructure interface loads are recovered using the mode acceleration method. The dis- 
placement DRM has two partitions. The first consists of the unit load static deflection 
vectors which are extracted from the static analysis database. The second, which provides 
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the dynamic contribution, is obtained from the normal modes database. A rigid format 
alter in the SOL 2 data recovery run assembles the full displacement DRM for responses 
specified through XYPLOT/XYPEAK requests. 

Displacement DRM = [rj) ^e/A“ ; a ] (8) 

The size of the DRMs is controlled by many factors. The number of rows in the acceleration 
and displacement DRMs will correspond to the number of response requests in the data 
recovery runs. Since the number of output requests is usually much smaller than the 
model size, data recovery operations using DRM procedures are much faster than standard 
methods. The number of columns in the acceleration DRM will correspond to the number 
of retained system modes. The number of columns of the displacement DRM will be equal 
to the sum of the retained elastic modes and the number of load application points. 

3.0 EZTRAN Execution: The solution of the generalized equations of motion and 
the recovery of physical responses are accomplished by EZTRAN using the NASTRAN 
generated data. The information provided to EZTRAN is shown in Figure 2. The gener- 
alized mass and stiffness matrices and the generalized unit forces, F u = (f> T I, are supplied 
by NASTRAN through a model file. The scaled forcing functions are supplied through 
a forcing function file. Specific instructions for an EZTRAN run including load cases to 
be analyzed, time step information, number of modes to be included, modal damping 
parameters, and type of initial conditions are entered by the user in an input file. The 
NASTRAN generated DRMs are supplied through a matrix file. A dictionary file provides 
identification for the response items in the DRMs. 

The modal equations of motion are uncoupled by virtue of linearity and proportional 
damping assumptions. 


\M\ {<?} + [CJ {<?} + \K] {(?} = {f„}{p(<)} (9) 

They are solved using a simple recursive algorithm [Rf. 4]. The solution is exact within the 
limits that the applied forces are assumed to vary linearly between integration steps. The 
method is unconditionally stable, regardless of integration step size. However, the step size 
must be sufficiently small so that linear interpolation accurately follows the applied force 
time histories. Initial conditions are either zero (undeformed structure) or can be auto- 
matically computed by EZTRAN, assuming that the system is in steady-state equilibrium 
with initial non-zero forces. For example, the Orbiter and payloads are initially deflected 
by gravity, wind loads, and restraining forces at the launch pad attach bolts, and are in 
steady state equilibrium. The deflections of elastic modes and acceleration of rigid body 
modes are computed from the initial modal forces by EZTRAN for such cases. 

The solution of modal differential equations yields the modal acceleration, Q. Physical 
responses are recovered using these solutions and the DRMs. 

f Physical 1 _ T Acceleration 1 r . \ 

\ Accelerations J ~ [ DRM J ( 10 ) 
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Physical Displacements and 
< Displacement dependent 
Responses 


Displacement 

DRM 



( 11 ) 


3.7 Post Processing: The results from EZTRAN include minimum/maximum sum- 
maries and time histories for the response items selected in the data recovery runs. The 
responses can be scaled by static and dynamic uncertainity factors to account for possible 
variations in the dynamic models or forcing functions. The results are written to format- 
ted files that can be read by other postprocessing programs to provide extrema reports, 
response history plots, shock response spectra, relative displacements, and other output. 

4. Example Problem: To illustrate the accuracy and efficiency of the alternate method, 
an example problem was analyzed. The problem consists of simple models of the Orbiter, 
Spacelab, Floor, Rack and a Box which were assembled into a solution system as shown in 
Figure 3. Modal reductions were performed at each stage of substructure assembly. The 
model was analyzed for a dynamic transient load case which had 45 load application points 
on the Orbiter substructure. 

The analysis was performed using three different approaches. The first analysis used 
the direct transient solution feature of NASTRAN to solve a full, unreduced model of the 
system (3971 DOF). Although this approach is not practical for most real world problems, 
it provides an accurate baseline solution without any modal truncation errors. The second 
analysis used the modal reduced system model (205 DOF) with modes up to 35 Hz being 
retained in the final solution system. The transient analysis was performed in the modal 
domain, and the physical responses were recovered using NASTRAN’s mode acceleration 
method. Finally, the transient analysis was performed on the same modal reduced system 
model using the alternate mode acceleration method. All three cases were undamped with 
zero initial conditions. The simulations were carried out for 0.5 seconds with an integration 
time step of 0.001 seconds. 

A comparison of the cost and accuracy of the three methods clearly demonstrates 
the merits of the alternate method. The axial forces in a CBAR element of the BOX 
substructure are shown in Figure 4. The alternate method produces results which are much 
closer to the baseline solution than the NASTRAN mode acceleration solution. Similar 
results were obtained for other displacement dependent responses. In addition to being 
accurate, the alternate method was also more efficient than the other methods, as shown in 
Figure 5. The computational advantages of the alternate method become more pronounced 
as the length of simulation increases. 

5. Conclusions: An accurate and efficient method for performing coupled transient loads 
analysis was presented and compared with the standard transient analysis capabilities of 
NASTRAN. The procedure uses data recovery matrices to reduce matrix size and compu- 
tation times. The mode acceleration method is incorporated in the DRM formulations to 
recover accurate displacements and displacement- dependent quantities like element stress- 
es, element forces and interface loads. The method is ideally suited for large, multi-level 
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substructured models. 
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ABSTRACT 

An algorithm for calculating acoustic intensities from a time-harmonic pressure field in an 
axisymmetric fluid region is presented. Acoustic pressures are computed in a mesh ofNASTRAN 
triangular finite elements of revolution (TRIAAX) using an analogy between the scalar wave 
equation and elasticity equations. Acoustic intensities are then calculated from pressures and 
pressure derivatives taken over the mesh of TRIAAX elements. Intensities are displayed as vectors 
indicating the directions and magnitudes of energy flow at all mesh points in the acoustic field. A 
prolate spheroidal shell is modeled with axisymmetric shell elements (CONEAX) and submerged 
in a fluid region of TRIAAX elements. The model is analyzed to illustrate the acoustic intensity 
method and the usefulness of energy flow paths in the understanding of the response of 
fluid-structure interaction problems. The structural-acoustic analogy used is summarized for 
completeness. This study uncovered a NASTRAN limitation involving numerical precision issues 
in the CONEAX stiffness calculation causing large errors in the system matrices for nearly 
cylindrical cones 


INTRODUCTION 

An acoustic intensity formulation for general, axisymmetric, fluid domains modeled with TRIAAX 
elements is presented. Numerical acoustic field solutions to fluid-structure interaction problems currently yield 
acoustic pressure fields, which may be used to locate high acoustic pressure concentrations. The motivation for 
calculatingand displaying acoustic intensities is to help visualize the energy flow paths which cause high pressure 
regions The energy (low fields can then help to identify dominant power paths which flow between structure 
and fluid, and therefore the important radiating parts of a structure. 

1 he general problem ol computing the interaction of an elastic structure with an acoustic fluid can be 
solved by combining a limte clement model of the structure with a fluid loading computed using boundary 
element [1-11], finite element [12-23], combined finite element/analytical [24-26], T-matrix [27-291 and 
approximate fluid loading [30-32] techniques. In the fluid finite element approach, the exterior fluid domain is 
modeled with finite elements truncated at a finite distance from the structure and terminated with an 
approximate radiation boundaiy condition to absorb outgoing waves. The principal computational trade-off 
between this approach and the boundary element approach is that the finite element approach yields large 
banded matrices, whereas the boundaiy element approach yields smaller, densely-populated matrices. This 
riade off sometimes favors the finite element approach for long, slender structures like ships which are 
naturally banded. In addition, only the fluid finite element approach has directly available an explicit fluid 
mesh which can be used tor graphical display of the wave motion through the fluid. Since a significant part of our 
interest involves the display of wave propagation through both structure and fluid, we therefore formulate this 
problem using the fluid finite element approach. The principal drawbacks to fluid finite element modeling are 
the need for an approximate radiation boundary condition at the outer fluid boundary, the requirements on mesh 
size and extent (sometimes leading to frequency-dependent fluid meshes [20]), and the difficulty of generating 
the fluid mesh. J b h 
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Direct frequency response NASTRAN [33] solutions for axisymmetric regions are described in general, 
and for domains that are defined using both structural and fluid subregions. Structural regions are modeled using 
standard structural finite elements (CONEAX, TRIAAX, TRAPAX); fluid regions are modeled with triangular 
elements of revolution (TRIAAX) using an analogy relating the Helmholtz equation to the elasticity equations 
used for the structural elements. The analogy is described in detail by Everstine [34] and is summarized for 
completeness here. The modeling of fluid-structure interaction between fluid and structure domains is also 
defined, as well as the application of acoustic boundary conditions to fluid models. 

The acoustic intensity formulation includes the definition of the intensity quantity, the algorithm used 
to calculate fluid particle velocities using the pressure distribution in a general triangularized domain, and the 
calculation of the acoustic intensity vector quantity from pressures and velocities. The formulation has been 
implemented in the program AcINT (Acoustic INTensity), which functions as a post-processor to NASTRAN. 
An example, a submerged prolate-spheroidal shell with two sets of boundary conditions, is analyzed for a ring 
load. The resulting acoustic intensity fields are displayed for a given excitation frequency to illustrate the energy 
flow paths which result. The acoustic intensity vector plots show the utility of the method in identifying dominant 
power paths in fluid-structure interaction problems. 

STRUCTURAL-ACOUSTIC ANALOGY 

From an engineering point of view, it is convenient to be able to make use of existing general purpose 
finite element codes for analyzing structural-acoustic problems. Finite element codes are widely available, 
versatile, reliable, well supported, and an abundance of pre- and postprocessors may be used with them. Thus we 
summarize in this section an analogy [34] between the equations of elasticity and the wave equation of acoustics. 
This analogy allows the coupled structural acoustic problem to be solved with standard finite element codes like 
NASTRAN. 


The pressure p in an acoustic field satisfies the wave equation 



( 1 ) 


where V 2 is the Laplacian operator, p is the dynamic fluid pressure, c is the wave speed, and dots denote partial 
differentiation with respect to time. 


On the other hand, the x-component of the Navier equations of elasticity, which are the equations 
solved by structural analysis computer programs, is 


A + 2G 
G 


U yy “t" Uj2 


A + G . 

G (V ^ 


+ W ,xz) + 




( 2 ) 


where u, v, and w are the Cartesian components of displacement, X is a Lame elastic constant, G is the shear 
modulus, f, is the x-component of body force per unit volume (e.g., gravity), q is the mass density, and commas 


denote partial differentiation. 


A comparison of Eqs. 1 and 2 indicates that elastic finite elements can be used to model scalar pressure 
fields if we let u, the x-component of displacement, represent p, set v = w = 0 everywhere, f, = 0, and X G. 
For three-dimensional analysis, the engineering constants consistent with this last requirement are 


E e 


10 20 G e , 



( 3 ) 


where the element shear modulus G e can be selected arbitrarily. The subscript “e” has been added to these 
constants to emphasize that they are merely numbers assigned to the elements. 

A variety of boundary conditions may also be imposed. At a pressure-release boundary, p = 0 is 
enforced explicitly like other displacement boundary conditions. For gradient conditions, the pressure gradient 
dp/dn is enforced at a boundary point by applying a “force” to the unconstrained DOF at that point equal to 

G e Adp/dn , where A is the area assigned to the point and n is the outward normal from the fluid region. For 

example, the plane wave absorbing boundary condition 

i£ = _£ (4) 

dn c 
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is enforced by applying to each point on the outer fluid boundary a “force” given by - (G t A/c)p . Since this 
force ’ is proportional to the first time derivative of the fundamental solution variable p, this boundary Condition 
is imposed in the analogy by attaching to the fluid DOF a “dashpot” of constant G e A/c . The Neumann condition 

dp/dn = 0 is the natural boundary condition under this analogy. The next higher order local radiation 
boundary condition, the curved wave absorbing boundary condition [23,36] 
d£ = _P _ P 
3n c r ’ 

where r is the radius of the boundary, is enforced under the analogy by attaching in parallel both a “dashpot” and a 
“spring” between each boundary point and ground. 


At a fluid-structure interface (an accelerating boundary), momentum and 
require that 


i£ 

dn 


~Q u„ 


continuity considerations 


( 6 ) 


where n is the normal at the interface, q is the mass density of the fluid, and u n is the normal component of fluid 

particle acceleration. Under the analogy, this condition is enforced by applying to the fluid DOF at a 
fluid-structure interface a “force” given by - (G e pA)ii n . 


Tb summarize, the wave equation, Eq. 1, can be solved with elastic finite elements if the 
three dimensional region is modeled with 3-D solid finite elements having material properties given by Eq. 3, 
and only one of the three Cartesian components of displacement is retained to represent the scalar variable p. In 
Cartesian coordinates, any of the three components can be used. The solution of axisymmetric problems in 
cylindrical coordinates follows the same approach except that the z-component of displacement is the only one 
which can be used to represent p. J 


FINITE ELEMENT FORMULATION OF FLUID-STRUCTURE INTERACTION 

. Th ere are two fundamental fluid-structure interaction problems of interest in structural acoustics: 
acoustic radiation, in which a submerged elastic body is subjected to a mechanical excitation applied to the 
structure, and scattering, in which the structure is subjected to an external incident pressure loading For general 
time-dependent problems, the excitation is an arbitraiy function of time, whereas in the time-harmonic case of 
interest here, the excitation has a single circular frequency u . 


Although our specific interest here is the time-harmonic case, we summarize the theory [22 35] for the 
more general case, which includes an incident loading as well. The radiation problem will be covered as a special 
case. For scattering, we assume, without loss of generality, that the incident wave propagates in the negative z 
direction. The speed of such propagation is c, the speed of sound in the fluid. 


Within the fluid region, the total pressure p satisfies the wave equation, Eq. 1. Since the incident 
lree held pressure p, is known, it is convenient to decompose the total pressure p into the sum of incident and 
scattered pressures 

P = Pi + P., 

each of which satisfies the wave equation. 


We now formulate the problem for finite element solution. Consider an arbitrary, submerged 
three-dimensional elastic structure subjected to either internal time-dependent loads or an external 
time-dependent incident pressure. If the structure is modeled with finite elements, the resulting matrix equation 
ol motion for the structural degrees of freedom (DOF) is 

Mu + Bti + Ku = F - GAp, (g) 

where M, B, and K are the structural mass, viscous damping, and stiffness matrices (dimension s x s), respectively 
u is the displacement vector for all structural DOF (wet and dry) in terms of the coordinate systems selected by 
the user (s x r), F is the vector of applied mechanical forces applied to the structure (s x r), G is the rectangular 
transformation matrix of direction cosines to transform a vector of outward normal forces at the wet points to a 
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vector of forces at all points in the coordinate systems selected by the user (sxf), A is the diagonal area matrix foi 
the wet surface (f x f), p is the vector of total fluid pressures (incident + scattered) applied at the wet grid points (1 
x r), and dots denote differentiation with respect to time. The pressure p is assumed positive in compression. In 
the' above dimensions, s denotes the total number of independent structural DOF (wet and dry), 1 denotes the 
number of fluid DOF (the number of wet points), and r denotes the number of load cases. If first order finite 
elements are used for the surface discretization, surface areas, normals, and the transformation matrix O can be 
obtained from the calculation of the load vector resulting from an outwardly directed static unit pressure load on 
the structure’s wet surface. The matrix product G A can then be interpreted as the matrix which converts a vector 
of negative fluid pressures to structural loads in the global coordinate system. The last two equations can be 
combined to yield 

Mu + Bu + Ku + GAp s = F - GAp,. 


A finite element model of the fluid region (with scattered pressure p, as the unknown) results in a 
matrix equation of the form 

Qp. + Cp, + Hp s = F*», < 10) 

where p, is the vector of scattered fluid pressures at the grid points of the fluid region, Q and H are the fluid 
“inertia" and “stiffness” matrices (analogous to M and K for structures), C is the “damping” matrix arising from 
the radiation boundary condition (Eq. 4), and F® is the “loading” applied to fluid DOF due to the flu id -structure 
interface condition, Eq. 6. Using the analogy described in the preceding section, elastic finite elements can be 
used to model both structural and fluid regions. Material constants assigned to the elastic elements used to 
model the fluid are specified according to Eq. 3. In three dimensions, elastic solid elements are used (e g., 
isoparametric bricks (IHEXi) for general 3-D analysis or solids of revolution (TRIAAX, TRAPAX) for 
axisymmetric analysis). 


At the fluid-structure interface, Eqs. 6 and 7 can be combined to yield 

dp, 

an 


e(u„, - u n ), 


(ii) 


where n is the outward unit normal, and u,„ and u„ arc. respectively, the incident and total outward normal 
components of fluid particle acceleration at the interface. Thus, from the analogy, we impose the fluid-structure 
interface condition by applying a “load” to each interface fluid point given by 

F (p) = - pG e A(u ni - ti n ), ^ 

where the first minus sign is introduced since, in the coupled problem, we choose n as the outward normal from 

the structure into the fluid, making n an inward normal for the fluid region. The normal displacements u„ are 

related to the total displacements u by the same rectangular transformation matrix G used above. 

r>r (13) 

u„ = G u, v ’ 

where the superscript T denotes the matrix transpose. Eqs. 10, 12, and 13 can be combined to yield 

Qp, + Cp, + Hp, - pG e (GA) T u = -pGeAiin,. (I 4 ) 

Since the fluid-structure coupling terms in Eqs. 9 and 14 are nonsymmetric, we symmetrize the problem [21] by 
using a new fluid unknown q such that 


= | Psdt, 


Ps 


(15) 


If Eq. 14 is integrated in time, and the fluid element “shear modulus” G c is chosen as 



the overall matrix system describing the coupled problem can be written as 


(16) 
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( 17 ) 
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where v m (-- u ni ) is the oulward normal component of incident fluid particle velocity. 


I he new variahlc q is, except lor a multiplicative constant, the velocity potential <p , since 

P = -Q<t>- (1 8) 

Lq. 17 could also be recast in terms of (p rather than q as the fundamental fluid unknown, but no particular 
advantage would result. In fact, the use ol q rather than <p has the practical advantage that the fluid pressure can 
be recovered directly from the finite element program as the time derivative (velocity) of the unknown q. 

lo summarize, both structural and fluid regions arc modeled with finite elements. For the fluid region, 
the material constants assigned to the finite elements are 

n _ - 1020 „ -1 -1 

~ > <J t = — , v e = unspecified, p c = — -, , in , 

Q Q qc 2 ( 19 ) 

where E e , G e , v„ and Q t are the Young’s modulus, shear modulus, Poisson’s ratio, and mass density, 
respectively, assigned to the fluid finite elements. The properties q and c above are the actual density and sound 
speed for the fluid medium. Die radiation boundary condition used is the plane wave approximation, Eq. 4, which 
appears to be adequate if the outer fluid boundary is sufficiently far from the structure [20]. With this boundary 
condition, matrix C in Eq. 17 arises from dashpots applied at the outer fluid boundary with damping constant 
- A/ (pc) at each grid point to which the area A has been assigned. At the fluid-structure interface, matrix GAis 
entered using the areas (or areal direction cosines) assigned to each wet degree of freedom. (Recall that GA can 
be interpreted as the matrix which converts a vector of negative fluid pressures to structural loads in the global 
coordinate system.) 


For radiation problems, the right-hand side of Eq. 17 can be simplified further since (he incident 
pressure p, is zero, and we obtain 
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( 20 ) 


We note that the structural and fluid unknowns are not sequenced as perhaps implied by the partitioned 
form ol Eq. 20. I he coupling matrix GA is quite sparse and has nonzeros only for matrix rows associated with the 
structural DOF at the fluid-structure interface and columns associated with the coincident fluid points. Thus, the 
grid points should be sequenced for minimum matrix bandwidth or profile as if the structural and fluid meshes 
comprised a single large mesh. As a result, the structural and fluid grid points will, in general, be interspersed in 
their numbering, and the system matrices will be sparse and banded. 

ACOUSTIC INTENSITY CALCULATIONS 

I he procedure for solving for the acoustic intensity field in an axisymmetric fluid finite element model 
using NASTRAN [33] and the acoustic intensity post-processor AcINT is: 


• run NASTRAN on a dynamically loaded finite element model of structural elements 
and fluid elements using direct frequency response analysis, and generate resultant 
nodal pressures for the fluid region(s); 

• run AcINT using the output from NASTRAN to calculate nodal fluid velocities and 
acoustic intensities. 

The nodal pressures, actually the cosine coefficients of nodal pressures, are computed by NASTRAN 
in response to the AXISYMMETRIC = COSINE command in the case control part of the input data. Since only 
cosine coefficients are requested, the 2, 4, and 6 DOF are removed by NASTRAN for harmonic zero The 1 
and 5 DOF must be constrained by the user, leaving the 3 DOF to represent the scalar fluid velocity potential. 
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The velocity potential integrated over time is the fluid pressure. Therefore, m the analogy, noda pressures 
obtained by the case control command VELOCITY = ALL (or SETID). The resultant fluid nodal pressures are 
used as illustrated below to compute acoustic velocities and intensities. 


Acoustic Intensity 

Intensity is defined as the time averaged product of a pressure with the in-phase component of particle 
velocity. For time-harmonic analysis, where complex numbers are used, this calculation may be visualized as 
taking the dot product of the pressure and velocity phasors. Multiplying one complex number by the in-phase 
part of another complex number is the same operation as multiplying the first number by the complex conjugate 
of the other number and taking the real part of the result: 

( 21 ) 

I = R [pv ], 


where p is pressure, and v is the complex conjugate of velocity. (Often, a factor of 1/2 appears in intensity 
equations. However there is no factor 1/2 in the equations if the assumption is made that pressures and velocities 
are “effective” rms values rather than amplitudes. With this assumption, consistency is maintained, and there is 
no mixing of effective and peak quantities in this formulation.) 


Acoustic Velocities 

The derivation of acoustic velocities for axisymmetric problems is performed for the cosine : coefficients 

of the Fourier summation about the axis of rotational symmetry (z), where the r, z, and rotation about 6 DOF 

are active. Only the r and z variations of the scalar pressure field are used to calculate the acoustic velocity vector 

field. The particle velocity in a fluid domain is defined as: 

i _ i ,d p? dp- ( 22 ) 

v = — Vp = — (— i + — k) • 

(OQ (OQ dr dz 

where q is the fluid density, co is the circular frequency,/ is the square root of -1, and are the pressure 


derivatives in the r and z directions, and \ and k arc unit vectors in the rand z directions, respectively. A first 
order finite difference approximation of the pressure derivatives at the nodes in an individual 1 RIAAX element 

can be made as shown in Figure 1. The pressure differences between nodes are divided by the distances between 

nodes to approximate the first derivative of pressure in the direction of the two nodes. The approximate pressure 



Figure 1. Pressure Gradient Approximation. 

gradient equations are written for all nodes connected to the node for which the velocity vector is to be calculated. 
In the case above, all elements connected to node 1 must be found, and equations are written for each node 
connected to node 1. No duplicate equations are written when two elements share a common edge. An 
overdetermined system of equations for the pressure derivatives in the r and z directions is the result, with one 
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equation for each node connected to the primary node. The system of equations is 



where Sj, and s^ are the r and z components of the unit vector from the primary node to the connecting node 
i, and N is the number of connecting nodes. 

The derivatives are determined approximately using a least squares approach. The particle velocity 
vector is then solved for using Equation 22, and the acoustic intensity vector is given by Equation 21. The 
axisymmetric acoustic intensity field for a complete domain is found by repeating this procedure for each node 
in the domain. In this way, a complete energy flow solution can be derived from nodal pressures and element 
connectivities. 

EXAMPLES 

The axisymmetric model of a submerged, half prolate spheroidal shell is shown in Figure 2. The 
structure has a semi axes of 10 m and 5 m. The material is steel with a uniform thickness of 25 mm, a Young’s 
modulus (E) of 2.074 Ell N/m 2 , a density (p) of 7860 kg/m 3 , a Poisson’s ratio (v) of 0.3, and a material loss 

factor of 0.0. The frequency range of interest was 100 to 500 Hz. The problem was analyzed for harmonic zero, 
commonly referred to as the “breathing mode” of the domain, implying no variation in the solution field about 
the z axis. 



Figure 2. Axisymmetric Model of a Submerged, Half Prolate-spheroid. 
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As shown in the diagram, the structure is submerged in fluid, which was also modeled and interlaced 
with the structure. The fluid is seawater, with a density (p) ol 1025 kg/m 1 and a speed of sound (c) of 1500 m/s. 
Two sets of boundary conditions along the left, vertical fluid edge were applied, and are discussed in more detail 
below. 


A finite element (CONEAX) model of the problem is shown in Figure 3. The structure was modeled 
using axisymmetric conical shell elements. The plate thickness, steel material properties, and the frequency 
range determine an estimated average flexural wavelength [38] of 0.7 m at 500 Hz, which for a mesh requirement 
of about eight elements/wavelength translates to a structural mesh density of 0.087 m/element length. 

The outer fluid was modeled using TRIAAX elements. The seawater properties and frequency range 
determine an upper wavelength of 15.0 m and a lower wavelength of 3.0 m. The upper wavelength determines 
the location of the outer boundary of the fluid mesh, shown as A in Figure 2, as one wavelength from the 
structure, or 15.0 m. The lower wavelength determines the fluid mesh density, which for a minimum of 8 
elements/wavelength, specifies a fluid mesh size of 0.375 m/element length. The fluid and structural meshes 
are distinct at the fluid-structure boundary, but with coincident nodes. Fhese coincident nodes are coupled by 
area matrices which map fluid pressures to structural forces, as described earlier. Only one fluid DOF was 
assigned to each mesh point. 

The area matrices are input with DMIG cards, which apply area values to the damping, or B2PP matrix. 
Acurrent limitation of NASTRAN is the program performing nonsymmetric system matrix decompositions when 
MPC data is used with DMIG input. Despite the DMIG input being declared symmetric, when the MPC 
equations are used to obtain the BDD matrix, NASTRAN changes the matrix trailer to nonsymmetric. All matrix 
operations become nonsymmetric as a result, greatly increasing computer time. A sequence of AL'I’ER 
statements may be used to restore the trailer to its symmetric form. One such sequence (for the 1990 version 
of NASTRAN, Rigid Format 8) is shown below for BDD. 


ALTER 119$ 

DIAGONAL BDD/AVEC/*COLUMN */().$ 
ADD AVEC,/PVEC/(0.0,0.0)$ 

MERGE BDD,„,PVEC,/BDDSYM/-l//6$ 
EQUIV BDDSYM,BDD$ 


Make BDD trailer symmetric 

Vector of ones 

Vector of zeros (P-Vcc) 

Dummy merge 

BDD now symmetric 


This alter is inserted before the FRRD module. With the BDD symmetry flag restored, subsequent matrix 
operations will take advantage of the symmetry of the system, reducing the required computer time. This 
NASTRAN bug has been fixed by Gordon Chan of Unisys for the 1992 program release. 

Since the structural elements are about one-quarter the size of the desired fluid elements, some mesh 
transitioning in the fluid meshes was required from the mesh density of the structure (0.087 m/element length) 
to the mesh density of the fluid (0.375 m/elemcnt length). The fluid clement material properties were assigned 
according to Eq. 19. 


Boundary and Loading Conditions 

Two sets of boundary conditions were applied in the problem. In both cases, a pc impedance, or plane 
wave absorbing boundary condition was applied to the curved outer fluid boundary; and a point forcing function 
was applied at z = 6.0 in the positive r direction along the structure, as shown in Figure 2. I he structure was 
constrained in all degrees of freedom at the upper left end. The left vertical fluid boundary was modeled in two 
ways: ( 1 ) as a rigid wall, which reflects incident waves, and (2) as a pc impedance absorber. DMIG statements 

inputting -A/pc values for all boundary points simulated the pc plane wave absorber. 


Results 

The problem may be accurately solved for any set of frequencies from 100 to 500 Hz. Detailed results 
are presented here for the 100 Hz case for Cases 1 and 2. The numerical output is examined several ways: power 
input and power radiation calculations, structural displacements, acoustic pressure contour plots, and acoustic 
intensity vector plots. This set of output provides a nearly complete solution to the structural-acoustic problem. 
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Figure 3. Finite Element Mesh. 


Since hoi h materials are lossless, all power entering the structure must eventually be radiated into the 
llmd. 1 his means the power input by the forcing function must equal the power exiting the system through the 
plane-wave absorbing outer boundaries. Power input is defined as 

Fin = l' m v in , (24) 


where F, n is the radial complex input lorce (in this case unity), and v* n is the complex conjugate of the 
corresponding velocity at the force point. The power radiated through the absorbing boundary may be found 


by integrating the calculated acoustic intensity field normal to the boundary: 


I * n dA, 


where C denotes the boundary contour, n is the outward normal vector, and dA is the incremental area. This 
integration may be converted to a summation over all nodes on the boundary, where 

ri 

Fr.-ui ' n, A,. (26) 

i “ \ 

In this case, n, and A, arc calculated for each node on the boundary, and used with the intensity vector at each 

node to calculate t lie power leaving the system. The summation over all boundary nodes gives the total radiated 
power. 
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Radiated power is commonly represented by a spherical pressure field with a reference radius of one 
yard. If the outer boundary is in the far-field, the pressure at one yard is 

p,„,(dB) - 20lo g y^ + 120, <27) 

where 120 dB is added to calculate the pressure relative to 1 pPa. 

Tkblel lists the power quantities defined above, and the relative error calculated by subtracting the 
radiated power from the input power and normalizing to the input power. The error, about 2% for each case, 
is probably due to the curved outer boundary. The approximate boundary condition was of a plane-wave 
absorbing boundary. Since neither the boundary nor the radiated waves is perfectly plane, small reflections at 
the boundary can occur, causing the power balance to be slightly in error. The radiated power, and therefore 
the pressure at 1 yard, is higher for Case 2, and the error is either due to the additional left qc boundary s 
absorbing more power, or the boundary’s causing a small shift in the frequency response of the system. A full 
frequency sweep would be required to determine the effects of the left absorbing boundary. 


Case 

Power Input (W) 

Power Radiated (W) 

%Error 

Pressure@l yard 
(dB re: 1 pPa) 

i 

1.499E-8 

1.519E-8 

-1.33 

91.9 

2 

8.988E-8 

8.798E-8 

2.11 

99.5 


"Ihblel. Power Results for 100 Hz. 


A plot of the displaced shapes of the structure for both cases is shown in Figure 4. 1 he displacement 
field of the structure is complex, and a time dependent animation of the structural response is required to 
visualize fully the movement. The plots shown here are at a single phase angle in the displacement cycle (292. 5 ) 
and show that the change in boundary condition does not significantly alter the structural response. A sma 
phase shift has occurred, but the general shape is the same for both cases. The ring loading causes the 
discontinuity in the waveform at z = 6.0. This point is the source of the waves travelling to the left and right from 
the load. 

Acoustic intensity vector plots are superimposed on acoustic pressure contour plots in Figures 5-8. 
Figure 5 shows the entire field for Case 1, and Figure 6 is a close-up of the field near the structure. Figure 7 
is a plot representing Case 2, and Figure 8 is a closer view. A common pressure scale is used for all plots, with 
the letters on the contours corresponding to the pressure levels (dB re: 1 p.Pa). The vector lengths are 
proportional to the log of the intensity magnitudes. The log of intensity is used in the plots to overcome the 1/r 
decay in intensity magnitude with distance. Since this is an axisymmetric analysis at circumferential harmonic 
zero, the acoustic fields are constant for all angles about the z axis, and no net energy may pass through the lower 
z boundary. All intensity vectors along the lower boundary therefore have zero radial components. 

The plots for Case 1 (Figure 5 and Figure 6) show the highest levels of far-field pressures to be in the 
r and z directions, with values between contours I and J, or 70 to 75 dB. Near-field pressure peaks are indicated 
near the structure by D contours, or about 100 dB levels. Examining the far-field intensity vectors show a 
far-field condition (all acoustic energy directed outward) at the outer boundary, with the dominant energy flow 
paths in the r and z directions. The rigid boundary condition along the left wall is evidenced by the absence ot 
any outward z directed intensity component along it. The near-field intensity vectors shows an energy path that 
begins at the load point (z = 6.0) and branches to the right and left. At the right, or bottom of the structure, energy 
flows along the fluid-structure boundary before travelling away toward the far-field at the r = 0 boundary. To 
the left of the load point, energy re-enters and re-exits the structure twice before radiating outward at the z - 0 
boundary. The circulation of power to the left of the load point causes “false sources” to appear where the energy 
re-exits the structure. Examination of the entire intensity field identifies the load point as the original source 
of power though. 

For Case 2, the use of the pc boundary (Figure 7 and Figure 8) and at the left edge of the fluid causes 
the acoustic pressures in the r direction, or upper left of the fluid domain, to decrease significantly from the rigid 
boundary case, from 75 to 55 dB. The acoustic intensity field shows the reason for the decrease in pressure; the 
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vectors along the pc boundary now have a z-directed component, implying power exiting the system through the 
boundary. The overall radiated power is higher for Case 2, though, due to an increase in the radiated energy 
along the midsection of the outer boundary. The near-field intensity plot reveals that the circulating power flows 
between structure and fluid to the left of the load point have disappeared. The energy now flows along the 
fluid-structure boundary and radiates outward at the top of the structure. To the right of the load however, power 
now circulates. The dominant path is from the load point to the right; into and out of the structure; then along 
the structure until r = 0, where some power reenters the structure, and the rest radiates outward along the z axis. 


NASTRAN DIFFICULTIES ENCOUNTERED 

Two important limitations of NASTRAN became apparant during this study: the BDD damping matrix's 
being specified as unsymmetric regardless of the symmetric nature of fluid structure interaction and absorbing 
boundary data input by DMIG cards; and the formulation of stiffness coefficients for CONE AX elements. The 
BDD matrix trailer may be restored to symmetric using the ALTER statements outlined in the Example section. 
The difficulties with the CONEAX formulation are not easily fixed however. 

Stiffnesses for CONEAX elements are computed analytically by NASTRAN, and involve the inverse 
of Ar/Al .where Ar is the difference in radii and A1 is the total distance between the grids defining an element. 
For perfectly cylindrical shell elements with no variation in radii ( Ar/Al = 0), a different formulation is used 
to avoid a floating point error caused by a division by zero. However, no provision is made for small relative 
variation in radii ( Ar/Al == 0), and for a small range of elements the analytical computation is corrupted when 
computer precision limits are reached. Sometimes the error is so drastic that negative values are obtained for 
self term (diagonal) stiffnesses. The negative stiffnesses are reported to the user when NASTRAN checks the 
system matrices for singularities. However, sometimes the error may be drastic in the positive sense, i.e., 
stiffnesses orders of magnitude too large. No error would be reported to the user, and the final solution would 
be incorrect. 


CONEAX stiffness errors were encountered for the example described here at the upper left end of 
the structure, where the elements become nearly cylindrical. In this case, the stiffnesses of several of the near 
cylindrical elements were output and analyzed for accuracy. The two end elements were found to have large 
errors in stiffness. Io solve the problem, the radii of the element grid points were set equal, and the end of the 
structure was approximated as purely cylindrical. 

A Possible programming solution to the sensitivity of CONEAX stiffnesses to small relative differences 
in radii is to approximate nearly cylindrical regions as cylindrical. For example, if for a given element Ar/Al 

is below some specified tolerance e, the second grid radius is set equal to the first grid radius. The resulting model 
would be a stepwise approximation of the nearly cylindrical region. The chief problem is how to choose e. Studies 
would have to be performed on ranges of nearly cylindrical elements using different levels of computer precision 
to determine the accuracy limits on the analytical stiffness computation method. 


CONCLUSIONS 

The combination of structural displacement plots, pressure contours, and acoustic intensity vector fields 
all serve to reveal the complete state of a structural-acoustic problem. However, one component of the response 
is missing: the energy flow within the structure. The circulating energy along the structural-acoustic boundary 
indicated by the intensity plots show power flowing through the structure. A formulation similar to that for 
acoustic intensity can be performed for the structure; however more than one wave type must be considered. 
For the axisymmetric shells of revolution (CONEAX) used here, for example, both flexural (composed of both 
shear and moment waves) and longitudinal waveforms may transport energy through structures. Methods have 
been developed for general three-dimensional structural models of beams (BAR) and plates (QUAD4) [39], but 
have not yet been extended to axisymmetric problems. This additional analysis tool will help improve 
considerably the understanding of structural-acoustic, frequency response problems. 
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Undisplaced Shape 

Displaced Shape (rigid Boundary Condition) 
Displaced Shape (rho*c Boundary Condition) 



Figure 4. Structural Displacements for both Boundary Condition Cases (phase = 292.5 ) 






ligure 5. Acoustic Intensities and Pressures for Rigid Vertical Boundaty; 

I inos Denote Constant Pressure Contours; Vectors Denote Acoustic Intensities. 
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Figure 6. Close-up of Acoustic Intensities and Pressures for Rigid Vertical Boundary; 
Lines Denote Constant Pressure Contours; Vectors Denote Acoustic Intensities. 
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Figure 7. Acoustic Intensities and Pressures for Absorbing Vertical Boundary; 
Lines Denote Constant Pressure Contours; Vectors Denote Acoustic Intensities. 
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Figure 8. Close-up of Acoustic Intensities and Pressures for Absorbing Vertical boundary; 
Lines Denote Constant Pressure Contours; Vectors Denote Acoustic Intensities. 
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